Here we assess whether the faces used in the experiments conveyed the intended emotional expression.
# setwd(paste0(getwd(), "/ratings")) # set working subdir
# pilot1 must be treated separately because it used different faces
pilot1 <- read_csv("pilot1_ratings.csv") %>%
mutate(
exp = "pilot1",
participant = as.factor(subject),
question = as.factor(question),
emotion = factor(substr(pic_id, 1, 1), labels = c("anger", "neutral", "happiness")),
picID = as.factor(substr(pic_id, 2, 2))
) %>%
dplyr::select(exp, participant, question, picID, emotion, ratings = "rating")
exps <- c( # variable with experiment identifier
# participants x picID x question
rep(sub("_ratings.csv", "", list.files(pattern = ".csv")[1]), each = (30 * 8 * 3)),
rep(sub("_ratings.csv", "", list.files(pattern = ".csv")[2]), each = (30 * 8 * 3)),
rep(sub("_ratings.csv", "", list.files(pattern = ".csv")[4]), each = (16 * 8 * 3))
)
ratings <- list.files(pattern = ".csv")[c(1, 2, 4)] %>% # list of csv files in current directory (excluding pilot1)
map_df(~ read_csv(.)) %>%
mutate(
exp = as.factor(exps),
participant = as.factor(subject),
question = as.factor(question),
picID = as.factor(pic_id),
emotion = case_when(
picID == "1" | picID == "5" ~ "anger",
picID == "2" | picID == "10" ~ "neutral",
picID == "3" | picID == "15" ~ "happiness",
picID == "4" | picID == "20" ~ "disgust"
)
) %>%
dplyr::select(exp, participant, question, picID, emotion, ratings = "rating") %>%
bind_rows(., pilot1) %>%
mutate(
variability = case_when(
exp != "pilot1" & picID == "1" |
exp != "pilot1" & picID == "2" |
exp != "pilot1" & picID == "3" |
exp != "pilot1" & picID == "4" ~ "low",
exp != "pilot1" & picID == "5" |
exp != "pilot1" & picID == "10" |
exp != "pilot1" & picID == "15" |
exp != "pilot1" & picID == "20" ~ "high",
exp == "pilot1" ~ "NA"
)
) %>%
unite(condition, c(emotion, variability), sep = "_", remove = FALSE) %>%
mutate(
exp = factor(exp,
levels = c("pilot1", "pilot2", "exp1", "exp2") # sort pilots first
),
condition = recode(factor(condition),
"anger_low" = "anger, low variability",
"neutral_low" = "neutral, low variability",
"happiness_low" = "happiness, low variability",
"disgust_low" = "disgust, low variability",
"anger_high" = "anger, high variability",
"neutral_high" = "neutral, high variability",
"happiness_high" = "happiness, high variability",
"disgust_high" = "disgust, high variability",
"anger_NA" = "anger",
"neutral_NA" = "neutral",
"happiness_NA" = "happiness"
)
) %>%
dplyr::select(-picID)
rm(pilot1)Face ratings were collected at the end of all experiments. In Pilot 1, we used 9 identities and 2 emotional expressions, whereas in the other experiments we used 2 identities and 3 expressions (disgust as a filler, never regularly presented during trials).
summary.ratings <- summarySEwithin(
data = ratings,
measurevar = "ratings",
betweenvars = "exp",
withinvars = c("condition", "question"),
idvar = "participant"
) %>%
mutate(
CI.95.low = ratings - ci,
CI.95.high = ratings + ci
) %>%
dplyr::select(exp, condition, question, ratings, CI.95.low, CI.95.high) %>%
mutate(condition = factor(condition,
levels = c(
"neutral", "anger", "happiness",
"neutral, low variability", "anger, low variability", "happiness, low variability", "disgust, low variability",
"neutral, high variability", "anger, high variability", "happiness, high variability", "disgust, high variability"
)
)) %>%
arrange(exp, question, condition)
kable(summary.ratings, digits = 2)| exp | condition | question | ratings | CI.95.low | CI.95.high |
|---|---|---|---|---|---|
| pilot1 | neutral | angry | 4.00 | 3.73 | 4.27 |
| pilot1 | anger | angry | 7.02 | 6.72 | 7.33 |
| pilot1 | happiness | angry | 1.13 | 0.95 | 1.30 |
| pilot1 | neutral | happy | 3.18 | 2.97 | 3.40 |
| pilot1 | anger | happy | 1.44 | 1.20 | 1.68 |
| pilot1 | happiness | happy | 6.86 | 6.62 | 7.11 |
| pilot2 | neutral, low variability | angry | 3.02 | 2.17 | 3.88 |
| pilot2 | anger, low variability | angry | 5.96 | 5.20 | 6.72 |
| pilot2 | happiness, low variability | angry | 1.34 | 1.04 | 1.63 |
| pilot2 | disgust, low variability | angry | 3.40 | 2.35 | 4.44 |
| pilot2 | neutral, high variability | angry | 4.40 | 3.32 | 5.47 |
| pilot2 | anger, high variability | angry | 8.65 | 8.30 | 9.00 |
| pilot2 | happiness, high variability | angry | 1.27 | 1.00 | 1.54 |
| pilot2 | disgust, high variability | angry | 4.27 | 2.88 | 5.66 |
| pilot2 | neutral, low variability | disgusted | 3.09 | 2.07 | 4.10 |
| pilot2 | anger, low variability | disgusted | 4.65 | 3.74 | 5.55 |
| pilot2 | happiness, low variability | disgusted | 1.59 | 1.24 | 1.93 |
| pilot2 | disgust, low variability | disgusted | 7.71 | 6.90 | 8.52 |
| pilot2 | neutral, high variability | disgusted | 3.90 | 2.86 | 4.93 |
| pilot2 | anger, high variability | disgusted | 3.59 | 2.52 | 4.65 |
| pilot2 | happiness, high variability | disgusted | 1.59 | 1.23 | 1.94 |
| pilot2 | disgust, high variability | disgusted | 7.71 | 6.47 | 8.95 |
| pilot2 | neutral, low variability | happy | 3.52 | 2.66 | 4.39 |
| pilot2 | anger, low variability | happy | 1.71 | 1.21 | 2.21 |
| pilot2 | happiness, low variability | happy | 7.46 | 6.71 | 8.21 |
| pilot2 | disgust, low variability | happy | 2.27 | 1.47 | 3.08 |
| pilot2 | neutral, high variability | happy | 2.77 | 2.06 | 3.49 |
| pilot2 | anger, high variability | happy | 1.27 | 0.92 | 1.63 |
| pilot2 | happiness, high variability | happy | 7.96 | 7.32 | 8.61 |
| pilot2 | disgust, high variability | happy | 1.46 | 0.96 | 1.97 |
| exp1 | neutral, low variability | angry | 2.29 | 1.87 | 2.70 |
| exp1 | anger, low variability | angry | 5.92 | 5.13 | 6.72 |
| exp1 | happiness, low variability | angry | 1.52 | 1.09 | 1.96 |
| exp1 | disgust, low variability | angry | 3.86 | 3.11 | 4.61 |
| exp1 | neutral, high variability | angry | 3.59 | 3.05 | 4.12 |
| exp1 | anger, high variability | angry | 8.46 | 7.86 | 9.05 |
| exp1 | happiness, high variability | angry | 1.39 | 1.15 | 1.63 |
| exp1 | disgust, high variability | angry | 4.79 | 3.82 | 5.76 |
| exp1 | neutral, low variability | disgusted | 2.79 | 2.36 | 3.22 |
| exp1 | anger, low variability | disgusted | 4.22 | 3.52 | 4.92 |
| exp1 | happiness, low variability | disgusted | 2.06 | 1.36 | 2.76 |
| exp1 | disgust, low variability | disgusted | 7.52 | 6.65 | 8.39 |
| exp1 | neutral, high variability | disgusted | 3.39 | 2.78 | 4.00 |
| exp1 | anger, high variability | disgusted | 3.62 | 2.86 | 4.38 |
| exp1 | happiness, high variability | disgusted | 2.42 | 1.61 | 3.24 |
| exp1 | disgust, high variability | disgusted | 7.76 | 6.85 | 8.66 |
| exp1 | neutral, low variability | happy | 2.86 | 2.40 | 3.31 |
| exp1 | anger, low variability | happy | 2.02 | 1.63 | 2.41 |
| exp1 | happiness, low variability | happy | 8.36 | 7.93 | 8.79 |
| exp1 | disgust, low variability | happy | 1.96 | 1.56 | 2.35 |
| exp1 | neutral, high variability | happy | 2.16 | 1.77 | 2.55 |
| exp1 | anger, high variability | happy | 1.72 | 1.27 | 2.17 |
| exp1 | happiness, high variability | happy | 8.29 | 7.88 | 8.70 |
| exp1 | disgust, high variability | happy | 1.59 | 1.33 | 1.85 |
| exp2 | neutral, low variability | angry | 2.60 | 2.15 | 3.06 |
| exp2 | anger, low variability | angry | 6.44 | 5.87 | 7.00 |
| exp2 | happiness, low variability | angry | 1.27 | 0.91 | 1.63 |
| exp2 | disgust, low variability | angry | 3.80 | 3.13 | 4.47 |
| exp2 | neutral, high variability | angry | 4.14 | 3.43 | 4.85 |
| exp2 | anger, high variability | angry | 8.50 | 7.97 | 9.03 |
| exp2 | happiness, high variability | angry | 1.64 | 1.10 | 2.17 |
| exp2 | disgust, high variability | angry | 3.94 | 3.09 | 4.78 |
| exp2 | neutral, low variability | disgusted | 2.40 | 1.86 | 2.95 |
| exp2 | anger, low variability | disgusted | 4.54 | 3.90 | 5.18 |
| exp2 | happiness, low variability | disgusted | 1.57 | 1.26 | 1.87 |
| exp2 | disgust, low variability | disgusted | 7.80 | 7.06 | 8.55 |
| exp2 | neutral, high variability | disgusted | 3.77 | 2.98 | 4.56 |
| exp2 | anger, high variability | disgusted | 3.14 | 2.39 | 3.88 |
| exp2 | happiness, high variability | disgusted | 1.87 | 1.23 | 2.51 |
| exp2 | disgust, high variability | disgusted | 8.07 | 7.39 | 8.74 |
| exp2 | neutral, low variability | happy | 3.70 | 3.14 | 4.27 |
| exp2 | anger, low variability | happy | 1.90 | 1.58 | 2.23 |
| exp2 | happiness, low variability | happy | 7.87 | 7.22 | 8.52 |
| exp2 | disgust, low variability | happy | 1.50 | 1.29 | 1.72 |
| exp2 | neutral, high variability | happy | 3.00 | 2.30 | 3.70 |
| exp2 | anger, high variability | happy | 1.50 | 1.05 | 1.96 |
| exp2 | happiness, high variability | happy | 8.24 | 7.81 | 8.66 |
| exp2 | disgust, high variability | happy | 1.37 | 1.15 | 1.58 |
ratings.plot <-
ggplot(
ratings,
aes(
x = question,
y = ratings,
color = emotion,
fill = emotion
)
) +
geom_pirate(
bars = FALSE,
cis = TRUE,
lines = TRUE, lines_params = list(color = "black"),
points = TRUE, points_params = list(shape = 16, size = 5, alpha = .2),
violins = TRUE, violins_params = list(size = 1),
show.legend = TRUE
) +
scale_fill_viridis_d(option = "cividis") +
scale_color_viridis_d(option = "cividis") +
scale_y_continuous(
limits = c(1, 9),
breaks = seq(1, 9, 1)
) +
geom_hline(
yintercept = seq(1, 9, 1),
linetype = "dotted",
colour = "#999999",
size = .8,
alpha = .5
) +
labs(
x = "question",
y = "ratings"
) +
facet_wrap(~exp, scales = "free") +
theme_EmoSSR
ratings.plot# ### save plots to file
# # all plots in one figure
# save_plot(paste0(getwd(), "/ratings.svg"),
# ratings.plot,
# base_aspect_ratio = 1,
# base_height = 11,
# base_width = 17
# )For the statistical analyses, we will fit Bayesian multilevel models using brms, separately for each question. The constant effect will be emotion. Intercepts and slopes will be allowed to vary as a function of participant and picID.
# subset data
pilot1.Qangry.ratings <- ratings %>%
dplyr::select(-c(variability, emotion)) %>%
filter(exp == "pilot1", question == "angry") %>%
droplevels() %>%
mutate(condition = factor(condition,
levels = c("anger", "neutral", "happiness") # sort anger first
))
pilot1.Qangry.model <- brm(ratings ~ condition # constant effects
+ (condition | participant), # varying effects
data = pilot1.Qangry.ratings, # subset data
family = cumulative("probit"), # likelihood function
prior = weak.priors, # priors
sample_prior = TRUE, # draw samples from the prior distribution
inits = "0", # initial parameter values in the chains
control = list(
adapt_delta = .99, # parameters of NUTS sampler
max_treedepth = 15
),
chains = num.chains, # number of chains
iter = num.iter, # number of iterations
warmup = num.warmup, # number of warm-up samples
thin = num.thin, # thinning rate
algorithm = "sampling", # type of sampling algorithm
cores = num.chains, # number of processor cores to use when running chains in parallel
seed = seed.smorfia # RNG seed
)
saveRDS(pilot1.Qangry.model, file = "pilot1_Qangry_model.rds", compress = "gzip") # save results# subset data
pilot1.Qangry.ratings <- ratings %>%
dplyr::select(-c(variability, emotion)) %>%
filter(exp == "pilot1", question == "angry") %>%
droplevels() %>%
mutate(condition = factor(condition,
levels = c("anger", "neutral", "happiness") # sort anger first
))
# load model
pilot1.Qangry.model <- readRDS("pilot1_Qangry_model.rds")
# summary posterior distributions & diagnostics
summary.pilot1.Qangry.model <-
describe_posterior(pilot1.Qangry.model,
centrality = "mean",
ci = 0.95,
ci_method = "hdi",
test = NULL,
diagnostic = "all"
) %>%
as_tibble() %>%
separate(Parameter, c("coef", "Parameter"), sep = "_") %>%
dplyr::select(-coef) %>%
mutate(Parameter = recode(factor(Parameter),
"Intercept.1." = "Intercept[1]",
"Intercept.2." = "Intercept[2]",
"Intercept.3." = "Intercept[3]",
"Intercept.4." = "Intercept[4]",
"Intercept.5." = "Intercept[5]",
"Intercept.6." = "Intercept[6]",
"Intercept.7." = "Intercept[7]",
"Intercept.8." = "Intercept[8]",
"conditionneutral" = "anger minus neutral",
"conditionhappiness" = "anger minus happiness"
)) %>%
dplyr::select(Parameter, Mean, HDI95_Low = CI_low, HDI95_High = CI_high, Estimated_Sample_Size = ESS, MonteCarlo_StdErr = MCSE)
kable(summary.pilot1.Qangry.model, digits = 2)| Parameter | Mean | HDI95_Low | HDI95_High | Estimated_Sample_Size | MonteCarlo_StdErr |
|---|---|---|---|---|---|
| Intercept[1] | -3.57 | -4.02 | -3.12 | 10162 | 0 |
| Intercept[2] | -2.90 | -3.32 | -2.47 | 10259 | 0 |
| Intercept[3] | -2.42 | -2.83 | -2.02 | 10291 | 0 |
| Intercept[4] | -2.00 | -2.39 | -1.60 | 10680 | 0 |
| Intercept[5] | -1.24 | -1.61 | -0.88 | 10587 | 0 |
| Intercept[6] | -0.69 | -1.04 | -0.34 | 10987 | 0 |
| Intercept[7] | -0.18 | -0.51 | 0.16 | 10937 | 0 |
| Intercept[8] | 0.46 | 0.11 | 0.79 | 10941 | 0 |
| anger minus neutral | -1.91 | -2.37 | -1.48 | 9933 | 0 |
| anger minus happiness | -4.05 | -4.68 | -3.41 | 10251 | 0 |
# MCMC chains
pilot1.Qangry.model %>%
spread_draws(`b_.*`, regex = TRUE) %>%
dplyr::select(
"Chain" = ".chain",
"Intercept[1]" = "b_Intercept[1]",
"Intercept[2]" = "b_Intercept[2]",
"Intercept[3]" = "b_Intercept[3]",
"Intercept[4]" = "b_Intercept[4]",
"Intercept[5]" = "b_Intercept[5]",
"Intercept[6]" = "b_Intercept[6]",
"Intercept[7]" = "b_Intercept[7]",
"Intercept[8]" = "b_Intercept[8]",
"anger minus neutral" = "b_conditionneutral",
"anger minus happiness" = "b_conditionhappiness"
) %>%
as.data.frame() %>%
mcmc_trace(facet_args = list(ncol = 2)) +
geom_vline(
xintercept = seq(1, 9, 1),
linetype = "dotted",
colour = "#999999",
size = .8,
alpha = .5
) +
ggtitle("Chain convergence") +
theme_EmoSSR +
theme(legend.position = "none")### observed vs. predicted data ###
pilot1.Qangry.ratings %>%
mutate(predicted = posterior_predict(pilot1.Qangry.model)[1, ]) %>% # take only first row of sampled ratings (as example)
dplyr::rename(observed = ratings) %>%
gather(key = data_type, value = ratings, observed:predicted) %>%
ggplot(aes(
x = condition,
y = ratings,
color = condition,
fill = condition
)) +
geom_pirate(
bars = FALSE,
cis = TRUE,
lines = TRUE, lines_params = list(color = "black"),
points = TRUE, points_params = list(shape = 16, size = 5, alpha = .2),
violins = TRUE, violins_params = list(size = 1),
show.legend = FALSE
) +
scale_fill_viridis_d(option = "cividis") +
scale_color_viridis_d(option = "cividis") +
scale_y_continuous(
limits = c(1, 9),
breaks = seq(1, 9, 1)
) +
geom_hline(
yintercept = seq(1, 9, 1),
linetype = "dotted",
colour = "#999999",
size = .8,
alpha = .5
) +
facet_wrap(~data_type, scales = "free") +
xlab("") +
ggtitle("pilot 1, \"anger\" question ") +
theme_EmoSSRpilot1.Qangry.posterior.test <- pilot1.Qangry.model %>%
spread_draws(`b_.*`, regex = TRUE) %>%
dplyr::select(
"anger minus neutral" = "b_conditionneutral",
"anger minus happiness" = "b_conditionhappiness"
) %>%
gather(key = "condition", value = "ratings") %>%
mutate(condition = as.factor(condition))
# summary
summary.pilot1.Qangry.posterior.test <- pilot1.Qangry.posterior.test %>%
group_by(condition) %>%
dplyr::summarize(
mean = mean(ratings),
hdi.95.low = hdi(ratings)[1],
hdi.95.high = hdi(ratings)[2],
evidence.ratio = ifelse(mean < 0,
length(which(ratings < 0)) / length(which(ratings > 0)),
length(which(ratings > 0)) / length(which(ratings < 0))
)
) %>%
ungroup()
kable(summary.pilot1.Qangry.posterior.test, digits = 2)| condition | mean | hdi.95.low | hdi.95.high | evidence.ratio |
|---|---|---|---|---|
| anger minus happiness | -4.05 | -4.68 | -3.41 | Inf |
| anger minus neutral | -1.91 | -2.37 | -1.48 | Inf |
par(mfrow = c(1, 2))
for (i in 1:length(levels(pilot1.Qangry.posterior.test$condition))) {
temp.data <- filter(pilot1.Qangry.posterior.test, condition == levels(condition)[i])
plotPost(
pull(temp.data, ratings),
credMass = 0.95,
compVal = 0,
ROPE = NULL,
xlab = "",
col = "#b3cde0",
showCurve = FALSE,
cex = 1,
main = levels(temp.data$condition)[i]
)
}# subset data
pilot1.Qhappy.ratings <- ratings %>%
dplyr::select(-c(variability, emotion)) %>%
filter(exp == "pilot1", question == "happy") %>%
droplevels() %>%
mutate(condition = factor(condition,
levels = c("happiness", "neutral", "anger") # sort happiness first
))
pilot1.Qhappy.model <- brm(ratings ~ condition
+ (condition | participant),
data = pilot1.Qhappy.ratings,
family = cumulative("probit"),
prior = weak.priors,
sample_prior = TRUE,
inits = "0",
control = list(
adapt_delta = .99,
max_treedepth = 15
),
chains = num.chains,
iter = num.iter,
warmup = num.warmup,
thin = num.thin,
algorithm = "sampling",
cores = num.chains,
seed = seed.smorfia
)
saveRDS(pilot1.Qhappy.model, file = "pilot1_Qhappy_model.rds", compress = "gzip")# subset data
pilot1.Qhappy.ratings <- ratings %>%
dplyr::select(-c(variability, emotion)) %>%
filter(exp == "pilot1", question == "happy") %>%
droplevels() %>%
mutate(condition = factor(condition,
levels = c("happiness", "neutral", "anger") # sort happiness first
))
# load model
pilot1.Qhappy.model <- readRDS("pilot1_Qhappy_model.rds")
# summary posterior distributions & diagnostics
summary.pilot1.Qhappy.model <-
describe_posterior(pilot1.Qhappy.model,
centrality = "mean",
ci = 0.95,
ci_method = "hdi",
test = NULL,
diagnostic = "all"
) %>%
as_tibble() %>%
separate(Parameter, c("coef", "Parameter"), sep = "_") %>%
dplyr::select(-coef) %>%
mutate(Parameter = recode(factor(Parameter),
"Intercept.1." = "Intercept[1]",
"Intercept.2." = "Intercept[2]",
"Intercept.3." = "Intercept[3]",
"Intercept.4." = "Intercept[4]",
"Intercept.5." = "Intercept[5]",
"Intercept.6." = "Intercept[6]",
"Intercept.7." = "Intercept[7]",
"Intercept.8." = "Intercept[8]",
"conditionneutral" = "happiness minus neutral",
"conditionanger" = "happiness minus anger"
)) %>%
dplyr::select(Parameter, Mean, HDI95_Low = CI_low, HDI95_High = CI_high, Estimated_Sample_Size = ESS, MonteCarlo_StdErr = MCSE)
kable(summary.pilot1.Qhappy.model, digits = 2)| Parameter | Mean | HDI95_Low | HDI95_High | Estimated_Sample_Size | MonteCarlo_StdErr |
|---|---|---|---|---|---|
| Intercept[1] | -3.83 | -4.31 | -3.34 | 7737 | 0 |
| Intercept[2] | -3.22 | -3.70 | -2.75 | 7787 | 0 |
| Intercept[3] | -2.63 | -3.09 | -2.18 | 7790 | 0 |
| Intercept[4] | -2.11 | -2.55 | -1.67 | 7885 | 0 |
| Intercept[5] | -1.12 | -1.52 | -0.73 | 8214 | 0 |
| Intercept[6] | -0.73 | -1.10 | -0.35 | 8914 | 0 |
| Intercept[7] | 0.01 | -0.36 | 0.37 | 9142 | 0 |
| Intercept[8] | 0.79 | 0.42 | 1.18 | 9953 | 0 |
| happiness minus neutral | -2.54 | -2.97 | -2.12 | 10829 | 0 |
| happiness minus anger | -3.96 | -4.49 | -3.40 | 9539 | 0 |
# MCMC chains
pilot1.Qhappy.model %>%
spread_draws(`b_.*`, regex = TRUE) %>%
dplyr::select("Chain" = ".chain",
"Intercept[1]" = "b_Intercept[1]",
"Intercept[2]" = "b_Intercept[2]",
"Intercept[3]" = "b_Intercept[3]",
"Intercept[4]" = "b_Intercept[4]",
"Intercept[5]" = "b_Intercept[5]",
"Intercept[6]" = "b_Intercept[6]",
"Intercept[7]" = "b_Intercept[7]",
"Intercept[8]" = "b_Intercept[8]",
"happiness minus neutral" = "b_conditionneutral",
"happiness minus anger" = "b_conditionanger"
) %>%
as.data.frame() %>%
mcmc_trace(facet_args = list(ncol = 2)) +
geom_vline(
xintercept = seq(1, 9, 1),
linetype = "dotted",
colour = "#999999",
size = .8,
alpha = .5
) +
ggtitle("Chain convergence") +
theme_EmoSSR +
theme(legend.position = "none")### observed vs. predicted data ###
pilot1.Qhappy.ratings %>%
mutate(predicted = posterior_predict(pilot1.Qhappy.model)[1, ]) %>% # take only first row of sampled ratings (as example)
dplyr::rename(observed = ratings) %>%
gather(key = data_type, value = ratings, observed:predicted) %>%
ggplot(aes(x = condition,
y = ratings,
color = condition,
fill = condition)) +
geom_pirate(
bars = FALSE,
cis = TRUE,
lines = TRUE, lines_params = list(color = "black"),
points = TRUE, points_params = list(shape = 16, size = 5, alpha = .2),
violins = TRUE, violins_params = list(size = 1),
show.legend = FALSE
) +
scale_fill_viridis_d(option = "cividis") +
scale_color_viridis_d(option = "cividis") +
scale_y_continuous(limits = c(1, 9),
breaks = seq(1, 9, 1)) +
geom_hline(yintercept = seq(1, 9, 1),
linetype = "dotted",
colour = "#999999",
size = .8,
alpha = .5) +
facet_wrap( ~ data_type, scales = "free") +
xlab("") +
ggtitle("pilot 1, \"happiness\" question ") +
theme_EmoSSRpilot1.Qhappy.posterior.test <- pilot1.Qhappy.model %>%
spread_draws(`b_.*`, regex = TRUE) %>%
dplyr::select("happiness minus neutral" = "b_conditionneutral",
"happiness minus anger" = "b_conditionanger"
) %>%
gather(key = "condition", value = "ratings") %>%
mutate(condition = as.factor(condition))
# summary
summary.pilot1.Qhappy.posterior.test <- pilot1.Qhappy.posterior.test %>%
group_by(condition) %>%
dplyr::summarize(
mean = mean(ratings),
hdi.95.low = hdi(ratings)[1],
hdi.95.high = hdi(ratings)[2],
evidence.ratio = ifelse(mean < 0,
length(which(ratings < 0)) / length(which(ratings > 0)),
length(which(ratings > 0)) / length(which(ratings < 0))
)
) %>%
ungroup()
kable(summary.pilot1.Qhappy.posterior.test, digits = 2)| condition | mean | hdi.95.low | hdi.95.high | evidence.ratio |
|---|---|---|---|---|
| happiness minus anger | -3.96 | -4.49 | -3.40 | Inf |
| happiness minus neutral | -2.54 | -2.97 | -2.12 | Inf |
par(mfrow = c(1, 2))
for(i in 1:length(levels(pilot1.Qhappy.posterior.test$condition))) {
temp.data <- filter(pilot1.Qhappy.posterior.test, condition == levels(condition)[i])
plotPost(
pull(temp.data, ratings),
credMass = 0.95,
compVal = 0,
ROPE = NULL,
xlab = "",
col = "#b3cde0",
showCurve = FALSE,
cex = 1,
main = levels(temp.data$condition)[i]
)
}# subset data
pilot2.Qangry.ratings <- ratings %>%
filter(exp == "pilot2", question == "angry") %>%
droplevels() %>%
mutate(emotion = factor(emotion,
levels = c("anger", "neutral", "happiness", "disgust") # sort anger first
))
pilot2.Qangry.model <- brm(ratings ~ emotion : variability
+ (condition | participant),
data = pilot2.Qangry.ratings,
family = cumulative("probit"),
prior = weak.priors,
sample_prior = TRUE,
inits = "0",
control = list(
adapt_delta = .99,
max_treedepth = 15
),
chains = num.chains,
iter = num.iter,
warmup = num.warmup,
thin = num.thin,
algorithm = "sampling",
cores = num.chains,
seed = seed.smorfia
)
saveRDS(pilot2.Qangry.model, file = "pilot2_Qangry_model.rds", compress = "gzip")# subset data
pilot2.Qangry.ratings <- ratings %>%
filter(exp == "pilot2", question == "angry") %>%
droplevels() %>%
mutate(emotion = factor(emotion,
levels = c("anger", "neutral", "happiness", "disgust") # sort anger first
))
# load model
pilot2.Qangry.model <- readRDS("pilot2_Qangry_model.rds")
# summary posterior distributions & diagnostics
summary.pilot2.Qangry.model <-
describe_posterior(pilot2.Qangry.model,
centrality = "mean",
ci = 0.95,
ci_method = "hdi",
test = NULL,
diagnostic = "all"
) %>%
as_tibble() %>%
separate(Parameter, c("coef", "Parameter"), sep = "_") %>%
dplyr::select(-coef) %>%
mutate(Parameter = recode(factor(Parameter),
"Intercept.1." = "Intercept[1]",
"Intercept.2." = "Intercept[2]",
"Intercept.3." = "Intercept[3]",
"Intercept.4." = "Intercept[4]",
"Intercept.5." = "Intercept[5]",
"Intercept.6." = "Intercept[6]",
"Intercept.7." = "Intercept[7]",
"Intercept.8." = "Intercept[8]",
"emotionanger.variabilityhigh" = "anger, high variability",
"emotionneutral.variabilityhigh" = "neutral, high variability",
"emotionhappiness.variabilityhigh" = "happiness, high variability",
"emotiondisgust.variabilityhigh" = "disgust, high variability",
"emotionanger.variabilitylow" = "anger, low variability",
"emotionneutral.variabilitylow" = "neutral, low variability",
"emotionhappiness.variabilitylow" = "happiness, low variability",
"emotiondisgust.variabilitylow" = "disgust, low variability"
)) %>%
dplyr::select(Parameter, Mean, HDI95_Low = CI_low, HDI95_High = CI_high, Estimated_Sample_Size = ESS, MonteCarlo_StdErr = MCSE)
kable(summary.pilot2.Qangry.model, digits = 2)| Parameter | Mean | HDI95_Low | HDI95_High | Estimated_Sample_Size | MonteCarlo_StdErr |
|---|---|---|---|---|---|
| Intercept[1] | -2.29 | -4.79 | 0.04 | 7398 | 0.01 |
| Intercept[2] | -0.65 | -3.08 | 1.59 | 8242 | 0.01 |
| Intercept[3] | 0.19 | -2.13 | 2.53 | 8451 | 0.01 |
| Intercept[4] | 0.83 | -1.47 | 3.23 | 8392 | 0.01 |
| Intercept[5] | 1.45 | -0.93 | 3.78 | 8145 | 0.01 |
| Intercept[6] | 2.79 | 0.27 | 5.22 | 7620 | 0.01 |
| Intercept[7] | 4.78 | 2.07 | 7.59 | 6950 | 0.02 |
| Intercept[8] | 7.24 | 4.25 | 10.47 | 7065 | 0.02 |
| anger, high variability | 7.44 | 4.34 | 10.38 | 6673 | 0.02 |
| neutral, high variability | 0.58 | -1.80 | 2.84 | 8635 | 0.01 |
| happiness, high variability | -4.55 | -7.37 | -1.82 | 8455 | 0.02 |
| disgust, high variability | 0.16 | -2.41 | 2.76 | 8806 | 0.01 |
| anger, low variability | 2.36 | 0.05 | 4.83 | 7834 | 0.01 |
| neutral, low variability | -0.99 | -3.32 | 1.24 | 7900 | 0.01 |
| happiness, low variability | -4.68 | -7.69 | -1.76 | 8966 | 0.02 |
| disgust, low variability | -0.34 | -2.79 | 2.04 | 9104 | 0.01 |
# MCMC chains
pilot2.Qangry.model %>%
spread_draws(`b_.*`, regex = TRUE) %>%
dplyr::select("Chain" = ".chain",
"Intercept[1]" = "b_Intercept[1]",
"Intercept[2]" = "b_Intercept[2]",
"Intercept[3]" = "b_Intercept[3]",
"Intercept[4]" = "b_Intercept[4]",
"Intercept[5]" = "b_Intercept[5]",
"Intercept[6]" = "b_Intercept[6]",
"Intercept[7]" = "b_Intercept[7]",
"Intercept[8]" = "b_Intercept[8]",
"anger, high variability" = "b_emotionanger:variabilityhigh",
"neutral, high variability" = "b_emotionneutral:variabilityhigh",
"happiness, high variability" = "b_emotionhappiness:variabilityhigh",
"disgust, high variability" = "b_emotiondisgust:variabilityhigh",
"anger, low variability" = "b_emotionanger:variabilitylow",
"neutral, low variability" = "b_emotionneutral:variabilitylow",
"happiness, low variability" = "b_emotionhappiness:variabilitylow",
"disgust, low variability" = "b_emotiondisgust:variabilitylow"
) %>%
as.data.frame() %>%
mcmc_trace(facet_args = list(ncol = 2)) +
geom_vline(xintercept = seq(1, 9, 1),
linetype = "dotted",
colour = "#999999",
size = .8,
alpha = .5
) +
ggtitle("Chain convergence") +
theme_EmoSSR +
theme(legend.position = "none")### observed vs. predicted data ###
pilot2.Qangry.ratings %>%
dplyr::select(exp, participant, question, emotion, variability, condition, ratings) %>% # reorder columns
mutate(predicted = posterior_predict(pilot2.Qangry.model)[1, ]) %>% # take only first row of sampled ratings (as example)
dplyr::rename(observed = ratings) %>%
gather(key = data_type, value = ratings, observed:predicted) %>%
ggplot(aes(x = condition,
y = ratings,
color = condition,
fill = condition)) +
geom_pirate(
bars = FALSE,
cis = TRUE,
lines = TRUE, lines_params = list(color = "black"),
points = TRUE, points_params = list(shape = 16, size = 5, alpha = .2),
violins = TRUE, violins_params = list(size = 1),
show.legend = FALSE
) +
scale_fill_viridis_d(option = "cividis") +
scale_color_viridis_d(option = "cividis") +
scale_y_continuous(limits = c(1, 9),
breaks = seq(1, 9, 1)) +
geom_hline(yintercept = seq(1, 9, 1),
linetype = "dotted",
colour = "#999999",
size = .8,
alpha = .5) +
facet_wrap( ~ data_type, scales = "free") +
xlab("") +
ggtitle("pilot 2, \"anger\" question ") +
theme_EmoSSRpilot2.Qangry.posterior.test <- pilot2.Qangry.model %>%
spread_draws(`b_.*`, regex = TRUE) %>%
dplyr::select(
"anger, high variability" = "b_emotionanger:variabilityhigh",
"neutral, high variability" = "b_emotionneutral:variabilityhigh",
"happiness, high variability" = "b_emotionhappiness:variabilityhigh",
"disgust, high variability" = "b_emotiondisgust:variabilityhigh",
"anger, low variability" = "b_emotionanger:variabilitylow",
"neutral, low variability" = "b_emotionneutral:variabilitylow",
"happiness, low variability" = "b_emotionhappiness:variabilitylow",
"disgust, low variability" = "b_emotiondisgust:variabilitylow"
) %>%
mutate(
### emotion (anger minus others)
## high variability
# anger minus neutral
diff.high.angerVSneutral = `anger, high variability` - `neutral, high variability`,
# anger minus happiness
diff.high.angerVShappiness = `anger, high variability` - `happiness, high variability`,
# anger minus disgust
diff.high.angerVSdisgust = `anger, high variability` - `disgust, high variability`,
## low variability
# anger minus neutral
diff.low.angerVSneutral = `anger, low variability` - `neutral, low variability`,
# anger minus happiness
diff.low.angerVShappiness = `anger, low variability` - `happiness, low variability`,
# anger minus disgust
diff.low.angerVSdisgust = `anger, low variability` - `disgust, low variability`,
### variability (high minus low)
# anger
diff.anger.highVSlow = `anger, high variability` - `anger, low variability`,
# neutral
diff.neutral.highVSlow = `neutral, high variability` - `neutral, low variability`,
# happiness
diff.happiness.highVSlow = `happiness, high variability` - `happiness, low variability`,
# disgust
diff.disgust.highVSlow = `disgust, high variability` - `disgust, low variability`
) %>%
gather(key = "condition", value = "ratings") %>%
filter(condition == c(
"diff.high.angerVSneutral",
"diff.high.angerVShappiness",
"diff.high.angerVSdisgust",
"diff.low.angerVSneutral",
"diff.low.angerVShappiness",
"diff.low.angerVSdisgust",
"diff.anger.highVSlow",
"diff.neutral.highVSlow",
"diff.happiness.highVSlow",
"diff.disgust.highVSlow"
)) %>%
mutate(condition = factor(condition,
levels = c(
"diff.high.angerVSneutral", "diff.high.angerVShappiness", "diff.high.angerVSdisgust",
"diff.low.angerVSneutral", "diff.low.angerVShappiness", "diff.low.angerVSdisgust",
"diff.anger.highVSlow", "diff.neutral.highVSlow", "diff.happiness.highVSlow", "diff.disgust.highVSlow"
)
))
# summary
summary.pilot2.Qangry.posterior.test <- pilot2.Qangry.posterior.test %>%
group_by(condition) %>%
dplyr::summarize(
mean = mean(ratings),
hdi.95.low = hdi(ratings)[1],
hdi.95.high = hdi(ratings)[2],
evidence.ratio = ifelse(mean < 0,
length(which(ratings < 0)) / length(which(ratings > 0)),
length(which(ratings > 0)) / length(which(ratings < 0))
)
) %>%
ungroup()
kable(summary.pilot2.Qangry.posterior.test, digits = 2)| condition | mean | hdi.95.low | hdi.95.high | evidence.ratio |
|---|---|---|---|---|
| diff.high.angerVSneutral | 6.86 | 4.55 | 9.27 | Inf |
| diff.high.angerVShappiness | 12.01 | 8.53 | 15.63 | Inf |
| diff.high.angerVSdisgust | 7.30 | 4.60 | 10.14 | Inf |
| diff.low.angerVSneutral | 3.36 | 1.69 | 4.98 | Inf |
| diff.low.angerVShappiness | 7.07 | 4.47 | 9.90 | Inf |
| diff.low.angerVSdisgust | 2.71 | 0.96 | 4.66 | 399.00 |
| diff.anger.highVSlow | 5.10 | 3.28 | 7.00 | Inf |
| diff.neutral.highVSlow | 1.56 | 0.24 | 2.80 | 113.29 |
| diff.happiness.highVSlow | 0.13 | -2.55 | 2.86 | 1.09 |
| diff.disgust.highVSlow | 0.52 | -1.35 | 2.27 | 2.77 |
### Anger vs. others ###
pilot2.Qangry.posterior.test.emoref <- filter(
pilot2.Qangry.posterior.test,
condition == c(
"diff.high.angerVSneutral", "diff.high.angerVShappiness", "diff.high.angerVSdisgust",
"diff.low.angerVSneutral", "diff.low.angerVShappiness", "diff.low.angerVSdisgust"
)
) %>%
droplevels()
par(mfrow = c(3, 2))
for(i in 1:length(levels(pilot2.Qangry.posterior.test.emoref$condition))) {
temp.data <- filter(pilot2.Qangry.posterior.test.emoref, condition == levels(condition)[i])
plotPost(
pull(temp.data, ratings),
credMass = 0.95,
compVal = 0,
ROPE = NULL,
xlab = "",
col = "#b3cde0",
showCurve = FALSE,
cex = 1,
main = levels(temp.data$condition)[i]
)
}### High vs. low variability ###
pilot2.Qangry.posterior.test.var <- filter(
pilot2.Qangry.posterior.test,
condition == c(
"diff.anger.highVSlow", "diff.neutral.highVSlow", "diff.happiness.highVSlow", "diff.disgust.highVSlow"
)
) %>%
droplevels()
par(mfrow = c(2, 2))
for(i in 1:length(levels(pilot2.Qangry.posterior.test.var$condition))) {
temp.data <- filter(pilot2.Qangry.posterior.test.var, condition == levels(condition)[i])
plotPost(
pull(temp.data, ratings),
credMass = 0.95,
compVal = 0,
ROPE = NULL,
xlab = "",
col = "#b3cde0",
showCurve = FALSE,
cex = 1,
main = levels(temp.data$condition)[i]
)
}# subset data
pilot2.Qhappy.ratings <- ratings %>%
filter(exp == "pilot2", question == "happy") %>%
droplevels() %>%
mutate(emotion = factor(emotion,
levels = c("happiness", "neutral", "anger", "disgust") # sort happiness first
))
pilot2.Qhappy.model <- brm(ratings ~ emotion : variability
+ (condition | participant),
data = pilot2.Qhappy.ratings,
family = cumulative("probit"),
prior = weak.priors,
sample_prior = TRUE,
inits = "0",
control = list(
adapt_delta = .99,
max_treedepth = 15
),
chains = num.chains,
iter = num.iter,
warmup = num.warmup,
thin = num.thin,
algorithm = "sampling",
cores = num.chains,
seed = seed.smorfia
)
saveRDS(pilot2.Qhappy.model, file = "pilot2_Qhappy_model.rds", compress = "gzip")# subset data
pilot2.Qhappy.ratings <- ratings %>%
filter(exp == "pilot2", question == "happy") %>%
droplevels() %>%
mutate(emotion = factor(emotion,
levels = c("happiness", "neutral", "anger", "disgust") # sort happiness first
))
# load model
pilot2.Qhappy.model <- readRDS("pilot2_Qhappy_model.rds")
# summary posterior distributions & diagnostics
summary.pilot2.Qhappy.model <-
describe_posterior(pilot2.Qhappy.model,
centrality = "mean",
ci = 0.95,
ci_method = "hdi",
test = NULL,
diagnostic = "all"
) %>%
as_tibble() %>%
separate(Parameter, c("coef", "Parameter"), sep = "_") %>%
dplyr::select(-coef) %>%
mutate(Parameter = recode(factor(Parameter),
"Intercept.1." = "Intercept[1]",
"Intercept.2." = "Intercept[2]",
"Intercept.3." = "Intercept[3]",
"Intercept.4." = "Intercept[4]",
"Intercept.5." = "Intercept[5]",
"Intercept.6." = "Intercept[6]",
"Intercept.7." = "Intercept[7]",
"Intercept.8." = "Intercept[8]",
"emotionanger.variabilityhigh" = "anger, high variability",
"emotionneutral.variabilityhigh" = "neutral, high variability",
"emotionhappiness.variabilityhigh" = "happiness, high variability",
"emotiondisgust.variabilityhigh" = "disgust, high variability",
"emotionanger.variabilitylow" = "anger, low variability",
"emotionneutral.variabilitylow" = "neutral, low variability",
"emotionhappiness.variabilitylow" = "happiness, low variability",
"emotiondisgust.variabilitylow" = "disgust, low variability"
)) %>%
dplyr::select(Parameter, Mean, HDI95_Low = CI_low, HDI95_High = CI_high, Estimated_Sample_Size = ESS, MonteCarlo_StdErr = MCSE)
kable(summary.pilot2.Qhappy.model, digits = 2)| Parameter | Mean | HDI95_Low | HDI95_High | Estimated_Sample_Size | MonteCarlo_StdErr |
|---|---|---|---|---|---|
| Intercept[1] | -1.89 | -4.15 | 0.44 | 8717 | 0.01 |
| Intercept[2] | -0.49 | -2.79 | 1.74 | 9417 | 0.01 |
| Intercept[3] | 0.69 | -1.60 | 3.03 | 9774 | 0.01 |
| Intercept[4] | 1.19 | -1.16 | 3.50 | 9573 | 0.01 |
| Intercept[5] | 2.39 | -0.06 | 4.86 | 9343 | 0.01 |
| Intercept[6] | 3.57 | 0.95 | 6.08 | 8820 | 0.01 |
| Intercept[7] | 5.45 | 2.53 | 8.31 | 7592 | 0.02 |
| Intercept[8] | 7.34 | 3.90 | 10.94 | 6094 | 0.02 |
| happiness, high variability | 6.00 | 3.17 | 9.15 | 7681 | 0.02 |
| neutral, high variability | -0.61 | -2.92 | 1.63 | 9815 | 0.01 |
| anger, high variability | -3.30 | -5.67 | -0.88 | 8862 | 0.01 |
| disgust, high variability | -3.66 | -6.67 | -0.89 | 9665 | 0.01 |
| happiness, low variability | 5.11 | 2.27 | 7.89 | 8506 | 0.02 |
| neutral, low variability | 0.20 | -2.16 | 2.43 | 9573 | 0.01 |
| anger, low variability | -2.16 | -4.53 | 0.15 | 9134 | 0.01 |
| disgust, low variability | -1.64 | -4.12 | 0.86 | 9796 | 0.01 |
# MCMC chains
pilot2.Qhappy.model %>%
spread_draws(`b_.*`, regex = TRUE) %>%
dplyr::select("Chain" = ".chain",
"Intercept[1]" = "b_Intercept[1]",
"Intercept[2]" = "b_Intercept[2]",
"Intercept[3]" = "b_Intercept[3]",
"Intercept[4]" = "b_Intercept[4]",
"Intercept[5]" = "b_Intercept[5]",
"Intercept[6]" = "b_Intercept[6]",
"Intercept[7]" = "b_Intercept[7]",
"Intercept[8]" = "b_Intercept[8]",
"anger, high variability" = "b_emotionanger:variabilityhigh",
"neutral, high variability" = "b_emotionneutral:variabilityhigh",
"happiness, high variability" = "b_emotionhappiness:variabilityhigh",
"disgust, high variability" = "b_emotiondisgust:variabilityhigh",
"anger, low variability" = "b_emotionanger:variabilitylow",
"neutral, low variability" = "b_emotionneutral:variabilitylow",
"happiness, low variability" = "b_emotionhappiness:variabilitylow",
"disgust, low variability" = "b_emotiondisgust:variabilitylow"
) %>%
as.data.frame() %>%
mcmc_trace(facet_args = list(ncol = 2)) +
geom_vline(xintercept = seq(1, 9, 1),
linetype = "dotted",
colour = "#999999",
size = .8,
alpha = .5
) +
ggtitle("Chain convergence") +
theme_EmoSSR +
theme(legend.position = "none")### observed vs. predicted data ###
pilot2.Qhappy.ratings %>%
dplyr::select(exp, participant, question, emotion, variability, condition, ratings) %>% # reorder columns
mutate(predicted = posterior_predict(pilot2.Qhappy.model)[1, ]) %>% # take only first row of sampled ratings (as example)
dplyr::rename(observed = ratings) %>%
gather(key = data_type, value = ratings, observed:predicted) %>%
ggplot(aes(x = condition,
y = ratings,
color = condition,
fill = condition)) +
geom_pirate(
bars = FALSE,
cis = TRUE,
lines = TRUE, lines_params = list(color = "black"),
points = TRUE, points_params = list(shape = 16, size = 5, alpha = .2),
violins = TRUE, violins_params = list(size = 1),
show.legend = FALSE
) +
scale_fill_viridis_d(option = "cividis") +
scale_color_viridis_d(option = "cividis") +
scale_y_continuous(limits = c(1, 9),
breaks = seq(1, 9, 1)) +
geom_hline(yintercept = seq(1, 9, 1),
linetype = "dotted",
colour = "#999999",
size = .8,
alpha = .5) +
facet_wrap( ~ data_type, scales = "free") +
xlab("") +
ggtitle("pilot 2, \"happiness\" question ") +
theme_EmoSSRpilot2.Qhappy.posterior.test <- pilot2.Qhappy.model %>%
spread_draws(`b_.*`, regex = TRUE) %>%
dplyr::select(
"anger, high variability" = "b_emotionanger:variabilityhigh",
"neutral, high variability" = "b_emotionneutral:variabilityhigh",
"happiness, high variability" = "b_emotionhappiness:variabilityhigh",
"disgust, high variability" = "b_emotiondisgust:variabilityhigh",
"anger, low variability" = "b_emotionanger:variabilitylow",
"neutral, low variability" = "b_emotionneutral:variabilitylow",
"happiness, low variability" = "b_emotionhappiness:variabilitylow",
"disgust, low variability" = "b_emotiondisgust:variabilitylow"
) %>%
mutate(
### emotion (happiness minus others)
## high variability
# happiness minus neutral
diff.high.happinessVSneutral = `happiness, high variability` - `neutral, high variability`,
# happiness minus anger
diff.high.happinessVSanger = `happiness, high variability` - `anger, high variability`,
# happiness minus disgust
diff.high.happinessVSdisgust = `happiness, high variability` - `disgust, high variability`,
## low variability
# happiness minus neutral
diff.low.happinessVSneutral = `happiness, low variability` - `neutral, low variability`,
# happiness minus anger
diff.low.happinessVSanger = `happiness, low variability` - `anger, low variability`,
# happiness minus disgust
diff.low.happinessVSdisgust = `happiness, low variability` - `disgust, low variability`,
### variability (high minus low)
# anger
diff.anger.highVSlow = `anger, high variability` - `anger, low variability`,
# neutral
diff.neutral.highVSlow = `neutral, high variability` - `neutral, low variability`,
# happiness
diff.happiness.highVSlow = `happiness, high variability` - `happiness, low variability`,
# disgust
diff.disgust.highVSlow = `disgust, high variability` - `disgust, low variability`
) %>%
gather(key = "condition", value = "ratings") %>%
filter(condition == c(
"diff.high.happinessVSneutral",
"diff.high.happinessVSanger",
"diff.high.happinessVSdisgust",
"diff.low.happinessVSneutral",
"diff.low.happinessVSanger",
"diff.low.happinessVSdisgust",
"diff.anger.highVSlow",
"diff.neutral.highVSlow",
"diff.happiness.highVSlow",
"diff.disgust.highVSlow"
)) %>%
mutate(condition = factor(condition,
levels = c(
"diff.high.happinessVSneutral", "diff.high.happinessVSanger", "diff.high.happinessVSdisgust",
"diff.low.happinessVSneutral", "diff.low.happinessVSanger", "diff.low.happinessVSdisgust",
"diff.anger.highVSlow", "diff.neutral.highVSlow", "diff.happiness.highVSlow", "diff.disgust.highVSlow"
)
))
# summary
summary.pilot2.Qhappy.posterior.test <- pilot2.Qhappy.posterior.test %>%
group_by(condition) %>%
dplyr::summarize(
mean = mean(ratings),
hdi.95.low = hdi(ratings)[1],
hdi.95.high = hdi(ratings)[2],
evidence.ratio = ifelse(mean < 0,
length(which(ratings < 0)) / length(which(ratings > 0)),
length(which(ratings > 0)) / length(which(ratings < 0))
)
) %>%
ungroup()
kable(summary.pilot2.Qhappy.posterior.test, digits = 2)| condition | mean | hdi.95.low | hdi.95.high | evidence.ratio |
|---|---|---|---|---|
| diff.high.happinessVSneutral | 6.59 | 4.09 | 9.26 | Inf |
| diff.high.happinessVSanger | 9.27 | 6.50 | 12.48 | Inf |
| diff.high.happinessVSdisgust | 9.65 | 6.46 | 13.26 | Inf |
| diff.low.happinessVSneutral | 4.95 | 3.08 | 7.23 | Inf |
| diff.low.happinessVSanger | 7.31 | 4.89 | 9.96 | Inf |
| diff.low.happinessVSdisgust | 6.75 | 4.30 | 9.53 | Inf |
| diff.anger.highVSlow | -1.13 | -2.55 | 0.12 | 19.00 |
| diff.neutral.highVSlow | -0.78 | -1.91 | 0.37 | 12.22 |
| diff.happiness.highVSlow | 0.89 | -0.80 | 2.72 | 6.77 |
| diff.disgust.highVSlow | -2.01 | -5.05 | 0.14 | 21.86 |
### Happiness vs. others ###
pilot2.Qhappy.posterior.test.emoref <- filter(
pilot2.Qhappy.posterior.test,
condition == c(
"diff.high.happinessVSneutral", "diff.high.happinessVSanger", "diff.high.happinessVSdisgust",
"diff.low.happinessVSneutral", "diff.low.happinessVSanger", "diff.low.happinessVSdisgust"
)
) %>%
droplevels()
par(mfrow = c(3, 2))
for(i in 1:length(levels(pilot2.Qhappy.posterior.test.emoref$condition))) {
temp.data <- filter(pilot2.Qhappy.posterior.test.emoref, condition == levels(condition)[i])
plotPost(
pull(temp.data, ratings),
credMass = 0.95,
compVal = 0,
ROPE = NULL,
xlab = "",
col = "#b3cde0",
showCurve = FALSE,
cex = 1,
main = levels(temp.data$condition)[i]
)
}### High vs. low variability ###
pilot2.Qhappy.posterior.test.var <- filter(
pilot2.Qhappy.posterior.test,
condition == c(
"diff.anger.highVSlow", "diff.neutral.highVSlow", "diff.happiness.highVSlow", "diff.disgust.highVSlow"
)
) %>%
droplevels()
par(mfrow = c(2, 2))
for(i in 1:length(levels(pilot2.Qhappy.posterior.test.var$condition))) {
temp.data <- filter(pilot2.Qhappy.posterior.test.var, condition == levels(condition)[i])
plotPost(
pull(temp.data, ratings),
credMass = 0.95,
compVal = 0,
ROPE = NULL,
xlab = "",
col = "#b3cde0",
showCurve = FALSE,
cex = 1,
main = levels(temp.data$condition)[i]
)
}# subset data
pilot2.Qdisgust.ratings <- ratings %>%
filter(exp == "pilot2", question == "disgusted") %>%
droplevels() %>%
mutate(emotion = factor(emotion,
levels = c("disgust", "neutral", "anger", "happiness") # sort disgust first
))
pilot2.Qdisgust.model <- brm(ratings ~ emotion : variability
+ (condition | participant),
data = pilot2.Qdisgust.ratings,
family = cumulative("probit"),
prior = weak.priors,
sample_prior = TRUE,
inits = "0",
control = list(
adapt_delta = .99,
max_treedepth = 15
),
chains = num.chains,
iter = num.iter,
warmup = num.warmup,
thin = num.thin,
algorithm = "sampling",
cores = num.chains,
seed = seed.smorfia
)
saveRDS(pilot2.Qdisgust.model, file = "pilot2_Qdisgust_model.rds", compress = "gzip")# subset data
pilot2.Qdisgust.ratings <- ratings %>%
filter(exp == "pilot2", question == "disgusted") %>%
droplevels() %>%
mutate(emotion = factor(emotion,
levels = c("disgust", "neutral", "anger", "happiness") # sort disgust first
))
# load model
pilot2.Qdisgust.model <- readRDS("pilot2_Qdisgust_model.rds")
# summary posterior distributions & diagnostics
summary.pilot2.Qdisgust.model <-
describe_posterior(pilot2.Qdisgust.model,
centrality = "mean",
ci = 0.95,
ci_method = "hdi",
test = NULL,
diagnostic = "all"
) %>%
as_tibble() %>%
separate(Parameter, c("coef", "Parameter"), sep = "_") %>%
dplyr::select(-coef) %>%
mutate(Parameter = recode(factor(Parameter),
"Intercept.1." = "Intercept[1]",
"Intercept.2." = "Intercept[2]",
"Intercept.3." = "Intercept[3]",
"Intercept.4." = "Intercept[4]",
"Intercept.5." = "Intercept[5]",
"Intercept.6." = "Intercept[6]",
"Intercept.7." = "Intercept[7]",
"Intercept.8." = "Intercept[8]",
"emotionanger.variabilityhigh" = "anger, high variability",
"emotionneutral.variabilityhigh" = "neutral, high variability",
"emotionhappiness.variabilityhigh" = "happiness, high variability",
"emotiondisgust.variabilityhigh" = "disgust, high variability",
"emotionanger.variabilitylow" = "anger, low variability",
"emotionneutral.variabilitylow" = "neutral, low variability",
"emotionhappiness.variabilitylow" = "happiness, low variability",
"emotiondisgust.variabilitylow" = "disgust, low variability"
)) %>%
dplyr::select(Parameter, Mean, HDI95_Low = CI_low, HDI95_High = CI_high, Estimated_Sample_Size = ESS, MonteCarlo_StdErr = MCSE)
kable(summary.pilot2.Qdisgust.model, digits = 2)| Parameter | Mean | HDI95_Low | HDI95_High | Estimated_Sample_Size | MonteCarlo_StdErr |
|---|---|---|---|---|---|
| Intercept[1] | -2.08 | -4.33 | 0.15 | 6373 | 0.01 |
| Intercept[2] | -0.95 | -3.22 | 1.18 | 6575 | 0.01 |
| Intercept[3] | -0.36 | -2.57 | 1.83 | 6661 | 0.01 |
| Intercept[4] | 0.14 | -2.11 | 2.29 | 6728 | 0.01 |
| Intercept[5] | 0.85 | -1.28 | 3.14 | 6702 | 0.01 |
| Intercept[6] | 1.41 | -0.85 | 3.58 | 6743 | 0.01 |
| Intercept[7] | 2.98 | 0.56 | 5.28 | 6567 | 0.01 |
| Intercept[8] | 4.92 | 2.05 | 7.64 | 6530 | 0.02 |
| disgust, high variability | 4.28 | 1.12 | 7.27 | 8458 | 0.02 |
| neutral, high variability | -0.37 | -2.67 | 1.77 | 7071 | 0.01 |
| anger, high variability | -0.61 | -2.81 | 1.55 | 6771 | 0.01 |
| happiness, high variability | -3.06 | -5.55 | -0.76 | 6634 | 0.02 |
| disgust, low variability | 3.47 | 0.92 | 5.99 | 6854 | 0.02 |
| neutral, low variability | -1.29 | -3.55 | 1.03 | 6857 | 0.01 |
| anger, low variability | 0.30 | -1.92 | 2.45 | 6670 | 0.01 |
| happiness, low variability | -2.74 | -5.18 | -0.44 | 6657 | 0.01 |
# MCMC chains
pilot2.Qdisgust.model %>%
spread_draws(`b_.*`, regex = TRUE) %>%
dplyr::select("Chain" = ".chain",
"Intercept[1]" = "b_Intercept[1]",
"Intercept[2]" = "b_Intercept[2]",
"Intercept[3]" = "b_Intercept[3]",
"Intercept[4]" = "b_Intercept[4]",
"Intercept[5]" = "b_Intercept[5]",
"Intercept[6]" = "b_Intercept[6]",
"Intercept[7]" = "b_Intercept[7]",
"Intercept[8]" = "b_Intercept[8]",
"anger, high variability" = "b_emotionanger:variabilityhigh",
"neutral, high variability" = "b_emotionneutral:variabilityhigh",
"happiness, high variability" = "b_emotionhappiness:variabilityhigh",
"disgust, high variability" = "b_emotiondisgust:variabilityhigh",
"anger, low variability" = "b_emotionanger:variabilitylow",
"neutral, low variability" = "b_emotionneutral:variabilitylow",
"happiness, low variability" = "b_emotionhappiness:variabilitylow",
"disgust, low variability" = "b_emotiondisgust:variabilitylow"
) %>%
as.data.frame() %>%
mcmc_trace(facet_args = list(ncol = 2)) +
geom_vline(xintercept = seq(1, 9, 1),
linetype = "dotted",
colour = "#999999",
size = .8,
alpha = .5
) +
ggtitle("Chain convergence") +
theme_EmoSSR +
theme(legend.position = "none")### observed vs. predicted data ###
pilot2.Qdisgust.ratings %>%
dplyr::select(exp, participant, question, emotion, variability, condition, ratings) %>% # reorder columns
mutate(predicted = posterior_predict(pilot2.Qdisgust.model)[1, ]) %>% # take only first row of sampled ratings (as example)
dplyr::rename(observed = ratings) %>%
gather(key = data_type, value = ratings, observed:predicted) %>%
ggplot(aes(x = condition,
y = ratings,
color = condition,
fill = condition)) +
geom_pirate(
bars = FALSE,
cis = TRUE,
lines = TRUE, lines_params = list(color = "black"),
points = TRUE, points_params = list(shape = 16, size = 5, alpha = .2),
violins = TRUE, violins_params = list(size = 1),
show.legend = FALSE
) +
scale_fill_viridis_d(option = "cividis") +
scale_color_viridis_d(option = "cividis") +
scale_y_continuous(limits = c(1, 9),
breaks = seq(1, 9, 1)) +
geom_hline(yintercept = seq(1, 9, 1),
linetype = "dotted",
colour = "#999999",
size = .8,
alpha = .5) +
facet_wrap( ~ data_type, scales = "free") +
xlab("") +
ggtitle("pilot 2, \"disgusted\" question ") +
theme_EmoSSRpilot2.Qdisgust.posterior.test <- pilot2.Qdisgust.model %>%
spread_draws(`b_.*`, regex = TRUE) %>%
dplyr::select(
"anger, high variability" = "b_emotionanger:variabilityhigh",
"neutral, high variability" = "b_emotionneutral:variabilityhigh",
"happiness, high variability" = "b_emotionhappiness:variabilityhigh",
"disgust, high variability" = "b_emotiondisgust:variabilityhigh",
"anger, low variability" = "b_emotionanger:variabilitylow",
"neutral, low variability" = "b_emotionneutral:variabilitylow",
"happiness, low variability" = "b_emotionhappiness:variabilitylow",
"disgust, low variability" = "b_emotiondisgust:variabilitylow"
) %>%
mutate(
### emotion (disgust minus others)
## high variability
# disgust minus neutral
diff.high.disgustVSneutral = `disgust, high variability` - `neutral, high variability`,
# disgust minus anger
diff.high.disgustVSanger = `disgust, high variability` - `anger, high variability`,
# disgust minus happiness
diff.high.disgustVShappiness = `disgust, high variability` - `happiness, high variability`,
## low variability
# disgust minus neutral
diff.low.disgustVSneutral = `disgust, low variability` - `neutral, low variability`,
# disgust minus anger
diff.low.disgustVSanger = `disgust, low variability` - `anger, low variability`,
# disgust minus happiness
diff.low.disgustVShappiness = `disgust, low variability` - `happiness, low variability`,
### variability (high minus low)
# anger
diff.anger.highVSlow = `anger, high variability` - `anger, low variability`,
# neutral
diff.neutral.highVSlow = `neutral, high variability` - `neutral, low variability`,
# happiness
diff.happiness.highVSlow = `happiness, high variability` - `happiness, low variability`,
# disgust
diff.disgust.highVSlow = `disgust, high variability` - `disgust, low variability`
) %>%
gather(key = "condition", value = "ratings") %>%
filter(condition == c(
"diff.high.disgustVSneutral",
"diff.high.disgustVSanger",
"diff.high.disgustVShappiness",
"diff.low.disgustVSneutral",
"diff.low.disgustVSanger",
"diff.low.disgustVShappiness",
"diff.anger.highVSlow",
"diff.neutral.highVSlow",
"diff.happiness.highVSlow",
"diff.disgust.highVSlow"
)) %>%
mutate(condition = factor(condition,
levels = c(
"diff.high.disgustVSneutral", "diff.high.disgustVSanger", "diff.high.disgustVShappiness",
"diff.low.disgustVSneutral", "diff.low.disgustVSanger", "diff.low.disgustVShappiness",
"diff.anger.highVSlow", "diff.neutral.highVSlow", "diff.happiness.highVSlow", "diff.disgust.highVSlow"
)
))
# summary
summary.pilot2.Qdisgust.posterior.test <- pilot2.Qdisgust.posterior.test %>%
group_by(condition) %>%
dplyr::summarize(
mean = mean(ratings),
hdi.95.low = hdi(ratings)[1],
hdi.95.high = hdi(ratings)[2],
evidence.ratio = ifelse(mean < 0,
length(which(ratings < 0)) / length(which(ratings > 0)),
length(which(ratings > 0)) / length(which(ratings < 0))
)
) %>%
ungroup()
kable(summary.pilot2.Qdisgust.posterior.test, digits = 2)| condition | mean | hdi.95.low | hdi.95.high | evidence.ratio |
|---|---|---|---|---|
| diff.high.disgustVSneutral | 4.66 | 2.20 | 7.56 | 799.00 |
| diff.high.disgustVSanger | 4.89 | 2.37 | 7.70 | Inf |
| diff.high.disgustVShappiness | 7.40 | 4.53 | 10.31 | Inf |
| diff.low.disgustVSneutral | 4.76 | 2.71 | 6.85 | Inf |
| diff.low.disgustVSanger | 3.19 | 1.65 | 4.81 | 1599.00 |
| diff.low.disgustVShappiness | 6.26 | 4.19 | 8.40 | Inf |
| diff.anger.highVSlow | -0.90 | -1.77 | -0.09 | 44.71 |
| diff.neutral.highVSlow | 0.91 | -0.30 | 2.16 | 17.60 |
| diff.happiness.highVSlow | -0.33 | -1.92 | 1.26 | 1.93 |
| diff.disgust.highVSlow | 0.81 | -1.67 | 3.48 | 2.90 |
### Disgust vs. others ###
pilot2.Qdisgust.posterior.test.emoref <- filter(
pilot2.Qdisgust.posterior.test,
condition == c(
"diff.high.disgustVSneutral", "diff.high.disgustVSanger", "diff.high.disgustVShappiness",
"diff.low.disgustVSneutral", "diff.low.disgustVSanger", "diff.low.disgustVShappiness"
)
) %>%
droplevels()
par(mfrow = c(3, 2))
for(i in 1:length(levels(pilot2.Qdisgust.posterior.test.emoref$condition))) {
temp.data <- filter(pilot2.Qdisgust.posterior.test.emoref, condition == levels(condition)[i])
plotPost(
pull(temp.data, ratings),
credMass = 0.95,
compVal = 0,
ROPE = NULL,
xlab = "",
col = "#b3cde0",
showCurve = FALSE,
cex = 1,
main = levels(temp.data$condition)[i]
)
}### High vs. low variability ###
pilot2.Qdisgust.posterior.test.var <- filter(
pilot2.Qdisgust.posterior.test,
condition == c(
"diff.anger.highVSlow", "diff.neutral.highVSlow", "diff.happiness.highVSlow", "diff.disgust.highVSlow"
)
) %>%
droplevels()
par(mfrow = c(2, 2))
for(i in 1:length(levels(pilot2.Qdisgust.posterior.test.var$condition))) {
temp.data <- filter(pilot2.Qdisgust.posterior.test.var, condition == levels(condition)[i])
plotPost(
pull(temp.data, ratings),
credMass = 0.95,
compVal = 0,
ROPE = NULL,
xlab = "",
col = "#b3cde0",
showCurve = FALSE,
cex = 1,
main = levels(temp.data$condition)[i]
)
}# subset data
exp1.Qangry.ratings <- ratings %>%
filter(exp == "exp1", question == "angry") %>%
droplevels() %>%
mutate(emotion = factor(emotion,
levels = c("anger", "neutral", "happiness", "disgust") # sort anger first
))
exp1.Qangry.model <- brm(ratings ~ emotion : variability
+ (condition | participant),
data = exp1.Qangry.ratings,
family = cumulative("probit"),
prior = weak.priors,
sample_prior = TRUE,
inits = "0",
control = list(
adapt_delta = .99,
max_treedepth = 15
),
chains = num.chains,
iter = num.iter,
warmup = num.warmup,
thin = num.thin,
algorithm = "sampling",
cores = num.chains,
seed = seed.smorfia
)
saveRDS(exp1.Qangry.model, file = "exp1_Qangry_model.rds", compress = "gzip")# subset data
exp1.Qangry.ratings <- ratings %>%
filter(exp == "exp1", question == "angry") %>%
droplevels() %>%
mutate(emotion = factor(emotion,
levels = c("anger", "neutral", "happiness", "disgust") # sort anger first
))
# load model
exp1.Qangry.model <- readRDS("exp1_Qangry_model.rds")
# summary posterior distributions & diagnostics
summary.exp1.Qangry.model <-
describe_posterior(exp1.Qangry.model,
centrality = "mean",
ci = 0.95,
ci_method = "hdi",
test = NULL,
diagnostic = "all"
) %>%
as_tibble() %>%
separate(Parameter, c("coef", "Parameter"), sep = "_") %>%
dplyr::select(-coef) %>%
mutate(Parameter = recode(factor(Parameter),
"Intercept.1." = "Intercept[1]",
"Intercept.2." = "Intercept[2]",
"Intercept.3." = "Intercept[3]",
"Intercept.4." = "Intercept[4]",
"Intercept.5." = "Intercept[5]",
"Intercept.6." = "Intercept[6]",
"Intercept.7." = "Intercept[7]",
"Intercept.8." = "Intercept[8]",
"emotionanger.variabilityhigh" = "anger, high variability",
"emotionneutral.variabilityhigh" = "neutral, high variability",
"emotionhappiness.variabilityhigh" = "happiness, high variability",
"emotiondisgust.variabilityhigh" = "disgust, high variability",
"emotionanger.variabilitylow" = "anger, low variability",
"emotionneutral.variabilitylow" = "neutral, low variability",
"emotionhappiness.variabilitylow" = "happiness, low variability",
"emotiondisgust.variabilitylow" = "disgust, low variability"
)) %>%
dplyr::select(Parameter, Mean, HDI95_Low = CI_low, HDI95_High = CI_high, Estimated_Sample_Size = ESS, MonteCarlo_StdErr = MCSE)
kable(summary.exp1.Qangry.model, digits = 2)| Parameter | Mean | HDI95_Low | HDI95_High | Estimated_Sample_Size | MonteCarlo_StdErr |
|---|---|---|---|---|---|
| Intercept[1] | -0.39 | -2.59 | 1.79 | 9138 | 0.01 |
| Intercept[2] | 0.10 | -2.09 | 2.26 | 9118 | 0.01 |
| Intercept[3] | 0.75 | -1.41 | 2.97 | 9129 | 0.01 |
| Intercept[4] | 1.29 | -0.90 | 3.50 | 8975 | 0.01 |
| Intercept[5] | 1.85 | -0.42 | 4.01 | 8797 | 0.01 |
| Intercept[6] | 2.36 | 0.11 | 4.58 | 8603 | 0.01 |
| Intercept[7] | 3.19 | 0.95 | 5.51 | 8169 | 0.01 |
| Intercept[8] | 4.20 | 1.86 | 6.58 | 7656 | 0.01 |
| anger, high variability | 4.27 | 1.96 | 6.61 | 7643 | 0.01 |
| neutral, high variability | 0.47 | -1.64 | 2.72 | 9162 | 0.01 |
| happiness, high variability | -3.06 | -5.98 | -0.32 | 10597 | 0.01 |
| disgust, high variability | 1.00 | -1.21 | 3.33 | 9176 | 0.01 |
| anger, low variability | 1.99 | -0.34 | 4.14 | 8742 | 0.01 |
| neutral, low variability | -0.71 | -3.03 | 1.40 | 9099 | 0.01 |
| happiness, low variability | -4.43 | -7.70 | -1.19 | 11415 | 0.02 |
| disgust, low variability | 0.48 | -1.69 | 2.70 | 9327 | 0.01 |
# MCMC chains
exp1.Qangry.model %>%
spread_draws(`b_.*`, regex = TRUE) %>%
dplyr::select("Chain" = ".chain",
"Intercept[1]" = "b_Intercept[1]",
"Intercept[2]" = "b_Intercept[2]",
"Intercept[3]" = "b_Intercept[3]",
"Intercept[4]" = "b_Intercept[4]",
"Intercept[5]" = "b_Intercept[5]",
"Intercept[6]" = "b_Intercept[6]",
"Intercept[7]" = "b_Intercept[7]",
"Intercept[8]" = "b_Intercept[8]",
"anger, high variability" = "b_emotionanger:variabilityhigh",
"neutral, high variability" = "b_emotionneutral:variabilityhigh",
"happiness, high variability" = "b_emotionhappiness:variabilityhigh",
"disgust, high variability" = "b_emotiondisgust:variabilityhigh",
"anger, low variability" = "b_emotionanger:variabilitylow",
"neutral, low variability" = "b_emotionneutral:variabilitylow",
"happiness, low variability" = "b_emotionhappiness:variabilitylow",
"disgust, low variability" = "b_emotiondisgust:variabilitylow"
) %>%
as.data.frame() %>%
mcmc_trace(facet_args = list(ncol = 2)) +
geom_vline(xintercept = seq(1, 9, 1),
linetype = "dotted",
colour = "#999999",
size = .8,
alpha = .5
) +
ggtitle("Chain convergence") +
theme_EmoSSR +
theme(legend.position = "none")### observed vs. predicted data ###
exp1.Qangry.ratings %>%
dplyr::select(exp, participant, question, emotion, variability, condition, ratings) %>% # reorder columns
mutate(predicted = posterior_predict(exp1.Qangry.model)[1, ]) %>% # take only first row of sampled ratings (as example)
dplyr::rename(observed = ratings) %>%
gather(key = data_type, value = ratings, observed:predicted) %>%
ggplot(aes(x = condition,
y = ratings,
color = condition,
fill = condition)) +
geom_pirate(
bars = FALSE,
cis = TRUE,
lines = TRUE, lines_params = list(color = "black"),
points = TRUE, points_params = list(shape = 16, size = 5, alpha = .2),
violins = TRUE, violins_params = list(size = 1),
show.legend = FALSE
) +
scale_fill_viridis_d(option = "cividis") +
scale_color_viridis_d(option = "cividis") +
scale_y_continuous(limits = c(1, 9),
breaks = seq(1, 9, 1)) +
geom_hline(yintercept = seq(1, 9, 1),
linetype = "dotted",
colour = "#999999",
size = .8,
alpha = .5) +
facet_wrap( ~ data_type, scales = "free") +
xlab("") +
ggtitle("pilot 2, \"anger\" question ") +
theme_EmoSSRexp1.Qangry.posterior.test <- exp1.Qangry.model %>%
spread_draws(`b_.*`, regex = TRUE) %>%
dplyr::select(
"anger, high variability" = "b_emotionanger:variabilityhigh",
"neutral, high variability" = "b_emotionneutral:variabilityhigh",
"happiness, high variability" = "b_emotionhappiness:variabilityhigh",
"disgust, high variability" = "b_emotiondisgust:variabilityhigh",
"anger, low variability" = "b_emotionanger:variabilitylow",
"neutral, low variability" = "b_emotionneutral:variabilitylow",
"happiness, low variability" = "b_emotionhappiness:variabilitylow",
"disgust, low variability" = "b_emotiondisgust:variabilitylow"
) %>%
mutate(
### emotion (anger minus others)
## high variability
# anger minus neutral
diff.high.angerVSneutral = `anger, high variability` - `neutral, high variability`,
# anger minus happiness
diff.high.angerVShappiness = `anger, high variability` - `happiness, high variability`,
# anger minus disgust
diff.high.angerVSdisgust = `anger, high variability` - `disgust, high variability`,
## low variability
# anger minus neutral
diff.low.angerVSneutral = `anger, low variability` - `neutral, low variability`,
# anger minus happiness
diff.low.angerVShappiness = `anger, low variability` - `happiness, low variability`,
# anger minus disgust
diff.low.angerVSdisgust = `anger, low variability` - `disgust, low variability`,
### variability (high minus low)
# anger
diff.anger.highVSlow = `anger, high variability` - `anger, low variability`,
# neutral
diff.neutral.highVSlow = `neutral, high variability` - `neutral, low variability`,
# happiness
diff.happiness.highVSlow = `happiness, high variability` - `happiness, low variability`,
# disgust
diff.disgust.highVSlow = `disgust, high variability` - `disgust, low variability`
) %>%
gather(key = "condition", value = "ratings") %>%
filter(condition == c(
"diff.high.angerVSneutral",
"diff.high.angerVShappiness",
"diff.high.angerVSdisgust",
"diff.low.angerVSneutral",
"diff.low.angerVShappiness",
"diff.low.angerVSdisgust",
"diff.anger.highVSlow",
"diff.neutral.highVSlow",
"diff.happiness.highVSlow",
"diff.disgust.highVSlow"
)) %>%
mutate(condition = factor(condition,
levels = c(
"diff.high.angerVSneutral", "diff.high.angerVShappiness", "diff.high.angerVSdisgust",
"diff.low.angerVSneutral", "diff.low.angerVShappiness", "diff.low.angerVSdisgust",
"diff.anger.highVSlow", "diff.neutral.highVSlow", "diff.happiness.highVSlow", "diff.disgust.highVSlow"
)
))
# summary
summary.exp1.Qangry.posterior.test <- exp1.Qangry.posterior.test %>%
group_by(condition) %>%
dplyr::summarize(
mean = mean(ratings),
hdi.95.low = hdi(ratings)[1],
hdi.95.high = hdi(ratings)[2],
evidence.ratio = ifelse(mean < 0,
length(which(ratings < 0)) / length(which(ratings > 0)),
length(which(ratings > 0)) / length(which(ratings < 0))
)
) %>%
ungroup()
kable(summary.exp1.Qangry.posterior.test, digits = 2)| condition | mean | hdi.95.low | hdi.95.high | evidence.ratio |
|---|---|---|---|---|
| diff.high.angerVSneutral | 3.79 | 2.86 | 4.79 | Inf |
| diff.high.angerVShappiness | 7.31 | 5.21 | 9.66 | Inf |
| diff.high.angerVSdisgust | 3.28 | 2.13 | 4.40 | Inf |
| diff.low.angerVSneutral | 2.73 | 1.73 | 3.69 | Inf |
| diff.low.angerVShappiness | 6.46 | 3.72 | 9.70 | Inf |
| diff.low.angerVSdisgust | 1.52 | 0.74 | 2.49 | Inf |
| diff.anger.highVSlow | 2.26 | 1.51 | 3.12 | Inf |
| diff.neutral.highVSlow | 1.18 | 0.39 | 1.92 | Inf |
| diff.happiness.highVSlow | 1.51 | -1.70 | 5.30 | 4.63 |
| diff.disgust.highVSlow | 0.53 | -0.39 | 1.46 | 6.11 |
### Anger vs. others ###
exp1.Qangry.posterior.test.emoref <- filter(
exp1.Qangry.posterior.test,
condition == c(
"diff.high.angerVSneutral", "diff.high.angerVShappiness", "diff.high.angerVSdisgust",
"diff.low.angerVSneutral", "diff.low.angerVShappiness", "diff.low.angerVSdisgust"
)
) %>%
droplevels()
par(mfrow = c(3, 2))
for(i in 1:length(levels(exp1.Qangry.posterior.test.emoref$condition))) {
temp.data <- filter(exp1.Qangry.posterior.test.emoref, condition == levels(condition)[i])
plotPost(
pull(temp.data, ratings),
credMass = 0.95,
compVal = 0,
ROPE = NULL,
xlab = "",
col = "#b3cde0",
showCurve = FALSE,
cex = 1,
main = levels(temp.data$condition)[i]
)
}### High vs. low variability ###
exp1.Qangry.posterior.test.var <- filter(
exp1.Qangry.posterior.test,
condition == c(
"diff.anger.highVSlow", "diff.neutral.highVSlow", "diff.happiness.highVSlow", "diff.disgust.highVSlow"
)
) %>%
droplevels()
par(mfrow = c(2, 2))
for(i in 1:length(levels(exp1.Qangry.posterior.test.var$condition))) {
temp.data <- filter(exp1.Qangry.posterior.test.var, condition == levels(condition)[i])
plotPost(
pull(temp.data, ratings),
credMass = 0.95,
compVal = 0,
ROPE = NULL,
xlab = "",
col = "#b3cde0",
showCurve = FALSE,
cex = 1,
main = levels(temp.data$condition)[i]
)
}# subset data
exp1.Qhappy.ratings <- ratings %>%
filter(exp == "exp1", question == "happy") %>%
droplevels() %>%
mutate(emotion = factor(emotion,
levels = c("happiness", "neutral", "anger", "disgust") # sort happiness first
))
exp1.Qhappy.model <- brm(ratings ~ emotion : variability
+ (condition | participant),
data = exp1.Qhappy.ratings,
family = cumulative("probit"),
prior = weak.priors,
sample_prior = TRUE,
inits = "0",
control = list(
adapt_delta = .99,
max_treedepth = 15
),
chains = num.chains,
iter = num.iter,
warmup = num.warmup,
thin = num.thin,
algorithm = "sampling",
cores = num.chains,
seed = seed.smorfia
)
saveRDS(exp1.Qhappy.model, file = "exp1_Qhappy_model.rds", compress = "gzip")# subset data
exp1.Qhappy.ratings <- ratings %>%
filter(exp == "exp1", question == "happy") %>%
droplevels() %>%
mutate(emotion = factor(emotion,
levels = c("happiness", "neutral", "anger", "disgust") # sort happiness first
))
# load model
exp1.Qhappy.model <- readRDS("exp1_Qhappy_model.rds")
# summary posterior distributions & diagnostics
summary.exp1.Qhappy.model <-
describe_posterior(exp1.Qhappy.model,
centrality = "mean",
ci = 0.95,
ci_method = "hdi",
test = NULL,
diagnostic = "all"
) %>%
as_tibble() %>%
separate(Parameter, c("coef", "Parameter"), sep = "_") %>%
dplyr::select(-coef) %>%
mutate(Parameter = recode(factor(Parameter),
"Intercept.1." = "Intercept[1]",
"Intercept.2." = "Intercept[2]",
"Intercept.3." = "Intercept[3]",
"Intercept.4." = "Intercept[4]",
"Intercept.5." = "Intercept[5]",
"Intercept.6." = "Intercept[6]",
"Intercept.7." = "Intercept[7]",
"Intercept.8." = "Intercept[8]",
"emotionanger.variabilityhigh" = "anger, high variability",
"emotionneutral.variabilityhigh" = "neutral, high variability",
"emotionhappiness.variabilityhigh" = "happiness, high variability",
"emotiondisgust.variabilityhigh" = "disgust, high variability",
"emotionanger.variabilitylow" = "anger, low variability",
"emotionneutral.variabilitylow" = "neutral, low variability",
"emotionhappiness.variabilitylow" = "happiness, low variability",
"emotiondisgust.variabilitylow" = "disgust, low variability"
)) %>%
dplyr::select(Parameter, Mean, HDI95_Low = CI_low, HDI95_High = CI_high, Estimated_Sample_Size = ESS, MonteCarlo_StdErr = MCSE)
kable(summary.exp1.Qhappy.model, digits = 2)| Parameter | Mean | HDI95_Low | HDI95_High | Estimated_Sample_Size | MonteCarlo_StdErr |
|---|---|---|---|---|---|
| Intercept[1] | -1.41 | -3.60 | 0.78 | 6491 | 0.01 |
| Intercept[2] | -0.41 | -2.58 | 1.79 | 6516 | 0.01 |
| Intercept[3] | 0.48 | -1.80 | 2.68 | 6386 | 0.01 |
| Intercept[4] | 0.84 | -1.35 | 3.14 | 6338 | 0.01 |
| Intercept[5] | 2.01 | -0.28 | 4.33 | 6408 | 0.01 |
| Intercept[6] | 2.46 | 0.17 | 4.83 | 6484 | 0.01 |
| Intercept[7] | 4.65 | 2.17 | 7.14 | 6622 | 0.02 |
| Intercept[8] | 6.58 | 3.74 | 9.46 | 6104 | 0.02 |
| happiness, high variability | 5.75 | 3.11 | 8.46 | 6326 | 0.02 |
| neutral, high variability | -1.63 | -3.92 | 0.55 | 6560 | 0.01 |
| anger, high variability | -2.41 | -4.63 | -0.18 | 6557 | 0.01 |
| disgust, high variability | -2.91 | -5.27 | -0.54 | 6647 | 0.01 |
| happiness, low variability | 5.79 | 3.15 | 8.49 | 6596 | 0.02 |
| neutral, low variability | -0.69 | -2.89 | 1.53 | 6549 | 0.01 |
| anger, low variability | -1.81 | -4.08 | 0.39 | 6539 | 0.01 |
| disgust, low variability | -2.25 | -4.62 | 0.08 | 6649 | 0.01 |
# MCMC chains
exp1.Qhappy.model %>%
spread_draws(`b_.*`, regex = TRUE) %>%
dplyr::select("Chain" = ".chain",
"Intercept[1]" = "b_Intercept[1]",
"Intercept[2]" = "b_Intercept[2]",
"Intercept[3]" = "b_Intercept[3]",
"Intercept[4]" = "b_Intercept[4]",
"Intercept[5]" = "b_Intercept[5]",
"Intercept[6]" = "b_Intercept[6]",
"Intercept[7]" = "b_Intercept[7]",
"Intercept[8]" = "b_Intercept[8]",
"anger, high variability" = "b_emotionanger:variabilityhigh",
"neutral, high variability" = "b_emotionneutral:variabilityhigh",
"happiness, high variability" = "b_emotionhappiness:variabilityhigh",
"disgust, high variability" = "b_emotiondisgust:variabilityhigh",
"anger, low variability" = "b_emotionanger:variabilitylow",
"neutral, low variability" = "b_emotionneutral:variabilitylow",
"happiness, low variability" = "b_emotionhappiness:variabilitylow",
"disgust, low variability" = "b_emotiondisgust:variabilitylow"
) %>%
as.data.frame() %>%
mcmc_trace(facet_args = list(ncol = 2)) +
geom_vline(xintercept = seq(1, 9, 1),
linetype = "dotted",
colour = "#999999",
size = .8,
alpha = .5
) +
ggtitle("Chain convergence") +
theme_EmoSSR +
theme(legend.position = "none")### observed vs. predicted data ###
exp1.Qhappy.ratings %>%
dplyr::select(exp, participant, question, emotion, variability, condition, ratings) %>% # reorder columns
mutate(predicted = posterior_predict(exp1.Qhappy.model)[1, ]) %>% # take only first row of sampled ratings (as example)
dplyr::rename(observed = ratings) %>%
gather(key = data_type, value = ratings, observed:predicted) %>%
ggplot(aes(x = condition,
y = ratings,
color = condition,
fill = condition)) +
geom_pirate(
bars = FALSE,
cis = TRUE,
lines = TRUE, lines_params = list(color = "black"),
points = TRUE, points_params = list(shape = 16, size = 5, alpha = .2),
violins = TRUE, violins_params = list(size = 1),
show.legend = FALSE
) +
scale_fill_viridis_d(option = "cividis") +
scale_color_viridis_d(option = "cividis") +
scale_y_continuous(limits = c(1, 9),
breaks = seq(1, 9, 1)) +
geom_hline(yintercept = seq(1, 9, 1),
linetype = "dotted",
colour = "#999999",
size = .8,
alpha = .5) +
facet_wrap( ~ data_type, scales = "free") +
xlab("") +
ggtitle("experiment 1, \"happiness\" question ") +
theme_EmoSSRexp1.Qhappy.posterior.test <- exp1.Qhappy.model %>%
spread_draws(`b_.*`, regex = TRUE) %>%
dplyr::select(
"anger, high variability" = "b_emotionanger:variabilityhigh",
"neutral, high variability" = "b_emotionneutral:variabilityhigh",
"happiness, high variability" = "b_emotionhappiness:variabilityhigh",
"disgust, high variability" = "b_emotiondisgust:variabilityhigh",
"anger, low variability" = "b_emotionanger:variabilitylow",
"neutral, low variability" = "b_emotionneutral:variabilitylow",
"happiness, low variability" = "b_emotionhappiness:variabilitylow",
"disgust, low variability" = "b_emotiondisgust:variabilitylow"
) %>%
mutate(
### emotion (happiness minus others)
## high variability
# happiness minus neutral
diff.high.happinessVSneutral = `happiness, high variability` - `neutral, high variability`,
# happiness minus anger
diff.high.happinessVSanger = `happiness, high variability` - `anger, high variability`,
# happiness minus disgust
diff.high.happinessVSdisgust = `happiness, high variability` - `disgust, high variability`,
## low variability
# happiness minus neutral
diff.low.happinessVSneutral = `happiness, low variability` - `neutral, low variability`,
# happiness minus anger
diff.low.happinessVSanger = `happiness, low variability` - `anger, low variability`,
# happiness minus disgust
diff.low.happinessVSdisgust = `happiness, low variability` - `disgust, low variability`,
### variability (high minus low)
# anger
diff.anger.highVSlow = `anger, high variability` - `anger, low variability`,
# neutral
diff.neutral.highVSlow = `neutral, high variability` - `neutral, low variability`,
# happiness
diff.happiness.highVSlow = `happiness, high variability` - `happiness, low variability`,
# disgust
diff.disgust.highVSlow = `disgust, high variability` - `disgust, low variability`
) %>%
gather(key = "condition", value = "ratings") %>%
filter(condition == c(
"diff.high.happinessVSneutral",
"diff.high.happinessVSanger",
"diff.high.happinessVSdisgust",
"diff.low.happinessVSneutral",
"diff.low.happinessVSanger",
"diff.low.happinessVSdisgust",
"diff.anger.highVSlow",
"diff.neutral.highVSlow",
"diff.happiness.highVSlow",
"diff.disgust.highVSlow"
)) %>%
mutate(condition = factor(condition,
levels = c(
"diff.high.happinessVSneutral", "diff.high.happinessVSanger", "diff.high.happinessVSdisgust",
"diff.low.happinessVSneutral", "diff.low.happinessVSanger", "diff.low.happinessVSdisgust",
"diff.anger.highVSlow", "diff.neutral.highVSlow", "diff.happiness.highVSlow", "diff.disgust.highVSlow"
)
))
# summary
summary.exp1.Qhappy.posterior.test <- exp1.Qhappy.posterior.test %>%
group_by(condition) %>%
dplyr::summarize(
mean = mean(ratings),
hdi.95.low = hdi(ratings)[1],
hdi.95.high = hdi(ratings)[2],
evidence.ratio = ifelse(mean < 0,
length(which(ratings < 0)) / length(which(ratings > 0)),
length(which(ratings > 0)) / length(which(ratings < 0))
)
) %>%
ungroup()
kable(summary.exp1.Qhappy.posterior.test, digits = 2)| condition | mean | hdi.95.low | hdi.95.high | evidence.ratio |
|---|---|---|---|---|
| diff.high.happinessVSneutral | 7.41 | 5.42 | 9.61 | Inf |
| diff.high.happinessVSanger | 8.15 | 5.89 | 10.20 | Inf |
| diff.high.happinessVSdisgust | 8.62 | 6.38 | 11.02 | Inf |
| diff.low.happinessVSneutral | 6.50 | 4.47 | 8.68 | Inf |
| diff.low.happinessVSanger | 7.60 | 5.67 | 9.93 | Inf |
| diff.low.happinessVSdisgust | 8.03 | 5.82 | 10.44 | Inf |
| diff.anger.highVSlow | -0.61 | -1.33 | 0.23 | 14.84 |
| diff.neutral.highVSlow | -0.93 | -1.79 | -0.12 | 176.78 |
| diff.happiness.highVSlow | -0.03 | -1.28 | 0.97 | 1.07 |
| diff.disgust.highVSlow | -0.68 | -2.08 | 0.79 | 6.17 |
### Happiness vs. others ###
exp1.Qhappy.posterior.test.emoref <- filter(
exp1.Qhappy.posterior.test,
condition == c(
"diff.high.happinessVSneutral", "diff.high.happinessVSanger", "diff.high.happinessVSdisgust",
"diff.low.happinessVSneutral", "diff.low.happinessVSanger", "diff.low.happinessVSdisgust"
)
) %>%
droplevels()
par(mfrow = c(3, 2))
for(i in 1:length(levels(exp1.Qhappy.posterior.test.emoref$condition))) {
temp.data <- filter(exp1.Qhappy.posterior.test.emoref, condition == levels(condition)[i])
plotPost(
pull(temp.data, ratings),
credMass = 0.95,
compVal = 0,
ROPE = NULL,
xlab = "",
col = "#b3cde0",
showCurve = FALSE,
cex = 1,
main = levels(temp.data$condition)[i]
)
}### High vs. low variability ###
exp1.Qhappy.posterior.test.var <- filter(
exp1.Qhappy.posterior.test,
condition == c(
"diff.anger.highVSlow", "diff.neutral.highVSlow", "diff.happiness.highVSlow", "diff.disgust.highVSlow"
)
) %>%
droplevels()
par(mfrow = c(2, 2))
for(i in 1:length(levels(exp1.Qhappy.posterior.test.var$condition))) {
temp.data <- filter(exp1.Qhappy.posterior.test.var, condition == levels(condition)[i])
plotPost(
pull(temp.data, ratings),
credMass = 0.95,
compVal = 0,
ROPE = NULL,
xlab = "",
col = "#b3cde0",
showCurve = FALSE,
cex = 1,
main = levels(temp.data$condition)[i]
)
}# subset data
exp1.Qdisgust.ratings <- ratings %>%
filter(exp == "exp1", question == "disgusted") %>%
droplevels() %>%
mutate(emotion = factor(emotion,
levels = c("disgust", "neutral", "anger", "happiness") # sort disgust first
))
exp1.Qdisgust.model <- brm(ratings ~ emotion : variability
+ (condition | participant),
data = exp1.Qdisgust.ratings,
family = cumulative("probit"),
prior = weak.priors,
sample_prior = TRUE,
inits = "0",
control = list(
adapt_delta = .99,
max_treedepth = 15
),
chains = num.chains,
iter = num.iter,
warmup = num.warmup,
thin = num.thin,
algorithm = "sampling",
cores = num.chains,
seed = seed.smorfia
)
saveRDS(exp1.Qdisgust.model, file = "exp1_Qdisgust_model.rds", compress = "gzip")# subset data
exp1.Qdisgust.ratings <- ratings %>%
filter(exp == "exp1", question == "disgusted") %>%
droplevels() %>%
mutate(emotion = factor(emotion,
levels = c("disgust", "neutral", "anger", "happiness") # sort disgust first
))
# load model
exp1.Qdisgust.model <- readRDS("exp1_Qdisgust_model.rds")
# summary posterior distributions & diagnostics
summary.exp1.Qdisgust.model <-
describe_posterior(exp1.Qdisgust.model,
centrality = "mean",
ci = 0.95,
ci_method = "hdi",
test = NULL,
diagnostic = "all"
) %>%
as_tibble() %>%
separate(Parameter, c("coef", "Parameter"), sep = "_") %>%
dplyr::select(-coef) %>%
mutate(Parameter = recode(factor(Parameter),
"Intercept.1." = "Intercept[1]",
"Intercept.2." = "Intercept[2]",
"Intercept.3." = "Intercept[3]",
"Intercept.4." = "Intercept[4]",
"Intercept.5." = "Intercept[5]",
"Intercept.6." = "Intercept[6]",
"Intercept.7." = "Intercept[7]",
"Intercept.8." = "Intercept[8]",
"emotionanger.variabilityhigh" = "anger, high variability",
"emotionneutral.variabilityhigh" = "neutral, high variability",
"emotionhappiness.variabilityhigh" = "happiness, high variability",
"emotiondisgust.variabilityhigh" = "disgust, high variability",
"emotionanger.variabilitylow" = "anger, low variability",
"emotionneutral.variabilitylow" = "neutral, low variability",
"emotionhappiness.variabilitylow" = "happiness, low variability",
"emotiondisgust.variabilitylow" = "disgust, low variability"
)) %>%
dplyr::select(Parameter, Mean, HDI95_Low = CI_low, HDI95_High = CI_high, Estimated_Sample_Size = ESS, MonteCarlo_StdErr = MCSE)
kable(summary.exp1.Qdisgust.model, digits = 2)| Parameter | Mean | HDI95_Low | HDI95_High | Estimated_Sample_Size | MonteCarlo_StdErr |
|---|---|---|---|---|---|
| Intercept[1] | -1.11 | -3.26 | 1.02 | 7706 | 0.01 |
| Intercept[2] | -0.41 | -2.57 | 1.68 | 7889 | 0.01 |
| Intercept[3] | 0.14 | -1.98 | 2.27 | 7925 | 0.01 |
| Intercept[4] | 0.70 | -1.43 | 2.82 | 7905 | 0.01 |
| Intercept[5] | 1.36 | -0.82 | 3.46 | 7833 | 0.01 |
| Intercept[6] | 1.90 | -0.22 | 4.11 | 7703 | 0.01 |
| Intercept[7] | 2.85 | 0.56 | 5.00 | 7245 | 0.01 |
| Intercept[8] | 4.42 | 2.08 | 6.81 | 6244 | 0.02 |
| disgust, high variability | 3.83 | 1.38 | 6.18 | 7161 | 0.01 |
| neutral, high variability | -0.47 | -2.60 | 1.68 | 8064 | 0.01 |
| anger, high variability | -0.17 | -2.33 | 1.89 | 8057 | 0.01 |
| happiness, high variability | -2.48 | -5.11 | 0.00 | 8998 | 0.01 |
| disgust, low variability | 3.16 | 0.82 | 5.51 | 7563 | 0.01 |
| neutral, low variability | -0.94 | -3.11 | 1.22 | 8027 | 0.01 |
| anger, low variability | 0.24 | -1.89 | 2.38 | 8179 | 0.01 |
| happiness, low variability | -3.22 | -5.83 | -0.59 | 9266 | 0.01 |
# MCMC chains
exp1.Qdisgust.model %>%
spread_draws(`b_.*`, regex = TRUE) %>%
dplyr::select("Chain" = ".chain",
"Intercept[1]" = "b_Intercept[1]",
"Intercept[2]" = "b_Intercept[2]",
"Intercept[3]" = "b_Intercept[3]",
"Intercept[4]" = "b_Intercept[4]",
"Intercept[5]" = "b_Intercept[5]",
"Intercept[6]" = "b_Intercept[6]",
"Intercept[7]" = "b_Intercept[7]",
"Intercept[8]" = "b_Intercept[8]",
"anger, high variability" = "b_emotionanger:variabilityhigh",
"neutral, high variability" = "b_emotionneutral:variabilityhigh",
"happiness, high variability" = "b_emotionhappiness:variabilityhigh",
"disgust, high variability" = "b_emotiondisgust:variabilityhigh",
"anger, low variability" = "b_emotionanger:variabilitylow",
"neutral, low variability" = "b_emotionneutral:variabilitylow",
"happiness, low variability" = "b_emotionhappiness:variabilitylow",
"disgust, low variability" = "b_emotiondisgust:variabilitylow"
) %>%
as.data.frame() %>%
mcmc_trace(facet_args = list(ncol = 2)) +
geom_vline(xintercept = seq(1, 9, 1),
linetype = "dotted",
colour = "#999999",
size = .8,
alpha = .5
) +
ggtitle("Chain convergence") +
theme_EmoSSR +
theme(legend.position = "none")### observed vs. predicted data ###
exp1.Qdisgust.ratings %>%
dplyr::select(exp, participant, question, emotion, variability, condition, ratings) %>% # reorder columns
mutate(predicted = posterior_predict(exp1.Qdisgust.model)[1, ]) %>% # take only first row of sampled ratings (as example)
dplyr::rename(observed = ratings) %>%
gather(key = data_type, value = ratings, observed:predicted) %>%
ggplot(aes(x = condition,
y = ratings,
color = condition,
fill = condition)) +
geom_pirate(
bars = FALSE,
cis = TRUE,
lines = TRUE, lines_params = list(color = "black"),
points = TRUE, points_params = list(shape = 16, size = 5, alpha = .2),
violins = TRUE, violins_params = list(size = 1),
show.legend = FALSE
) +
scale_fill_viridis_d(option = "cividis") +
scale_color_viridis_d(option = "cividis") +
scale_y_continuous(limits = c(1, 9),
breaks = seq(1, 9, 1)) +
geom_hline(yintercept = seq(1, 9, 1),
linetype = "dotted",
colour = "#999999",
size = .8,
alpha = .5) +
facet_wrap( ~ data_type, scales = "free") +
xlab("") +
ggtitle("experiment 1, \"disgusted\" question ") +
theme_EmoSSRexp1.Qdisgust.posterior.test <- exp1.Qdisgust.model %>%
spread_draws(`b_.*`, regex = TRUE) %>%
dplyr::select(
"anger, high variability" = "b_emotionanger:variabilityhigh",
"neutral, high variability" = "b_emotionneutral:variabilityhigh",
"happiness, high variability" = "b_emotionhappiness:variabilityhigh",
"disgust, high variability" = "b_emotiondisgust:variabilityhigh",
"anger, low variability" = "b_emotionanger:variabilitylow",
"neutral, low variability" = "b_emotionneutral:variabilitylow",
"happiness, low variability" = "b_emotionhappiness:variabilitylow",
"disgust, low variability" = "b_emotiondisgust:variabilitylow"
) %>%
mutate(
### emotion (disgust minus others)
## high variability
# disgust minus neutral
diff.high.disgustVSneutral = `disgust, high variability` - `neutral, high variability`,
# disgust minus anger
diff.high.disgustVSanger = `disgust, high variability` - `anger, high variability`,
# disgust minus happiness
diff.high.disgustVShappiness = `disgust, high variability` - `happiness, high variability`,
## low variability
# disgust minus neutral
diff.low.disgustVSneutral = `disgust, low variability` - `neutral, low variability`,
# disgust minus anger
diff.low.disgustVSanger = `disgust, low variability` - `anger, low variability`,
# disgust minus happiness
diff.low.disgustVShappiness = `disgust, low variability` - `happiness, low variability`,
### variability (high minus low)
# anger
diff.anger.highVSlow = `anger, high variability` - `anger, low variability`,
# neutral
diff.neutral.highVSlow = `neutral, high variability` - `neutral, low variability`,
# happiness
diff.happiness.highVSlow = `happiness, high variability` - `happiness, low variability`,
# disgust
diff.disgust.highVSlow = `disgust, high variability` - `disgust, low variability`
) %>%
gather(key = "condition", value = "ratings") %>%
filter(condition == c(
"diff.high.disgustVSneutral",
"diff.high.disgustVSanger",
"diff.high.disgustVShappiness",
"diff.low.disgustVSneutral",
"diff.low.disgustVSanger",
"diff.low.disgustVShappiness",
"diff.anger.highVSlow",
"diff.neutral.highVSlow",
"diff.happiness.highVSlow",
"diff.disgust.highVSlow"
)) %>%
mutate(condition = factor(condition,
levels = c(
"diff.high.disgustVSneutral", "diff.high.disgustVSanger", "diff.high.disgustVShappiness",
"diff.low.disgustVSneutral", "diff.low.disgustVSanger", "diff.low.disgustVShappiness",
"diff.anger.highVSlow", "diff.neutral.highVSlow", "diff.happiness.highVSlow", "diff.disgust.highVSlow"
)
))
# summary
summary.exp1.Qdisgust.posterior.test <- exp1.Qdisgust.posterior.test %>%
group_by(condition) %>%
dplyr::summarize(
mean = mean(ratings),
hdi.95.low = hdi(ratings)[1],
hdi.95.high = hdi(ratings)[2],
evidence.ratio = ifelse(mean < 0,
length(which(ratings < 0)) / length(which(ratings > 0)),
length(which(ratings > 0)) / length(which(ratings < 0))
)
) %>%
ungroup()
kable(summary.exp1.Qdisgust.posterior.test, digits = 2)| condition | mean | hdi.95.low | hdi.95.high | evidence.ratio |
|---|---|---|---|---|
| diff.high.disgustVSneutral | 4.31 | 2.94 | 5.98 | Inf |
| diff.high.disgustVSanger | 3.99 | 2.73 | 5.45 | Inf |
| diff.high.disgustVShappiness | 6.31 | 3.97 | 8.79 | Inf |
| diff.low.disgustVSneutral | 4.11 | 2.71 | 5.57 | Inf |
| diff.low.disgustVSanger | 2.93 | 1.66 | 4.17 | Inf |
| diff.low.disgustVShappiness | 6.37 | 4.15 | 8.85 | Inf |
| diff.anger.highVSlow | -0.41 | -1.07 | 0.25 | 7.60 |
| diff.neutral.highVSlow | 0.48 | -0.34 | 1.27 | 7.94 |
| diff.happiness.highVSlow | 0.73 | -1.46 | 2.74 | 3.32 |
| diff.disgust.highVSlow | 0.68 | -0.37 | 1.96 | 8.14 |
### Disgust vs. others ###
exp1.Qdisgust.posterior.test.emoref <- filter(
exp1.Qdisgust.posterior.test,
condition == c(
"diff.high.disgustVSneutral", "diff.high.disgustVSanger", "diff.high.disgustVShappiness",
"diff.low.disgustVSneutral", "diff.low.disgustVSanger", "diff.low.disgustVShappiness"
)
) %>%
droplevels()
par(mfrow = c(3, 2))
for(i in 1:length(levels(exp1.Qdisgust.posterior.test.emoref$condition))) {
temp.data <- filter(exp1.Qdisgust.posterior.test.emoref, condition == levels(condition)[i])
plotPost(
pull(temp.data, ratings),
credMass = 0.95,
compVal = 0,
ROPE = NULL,
xlab = "",
col = "#b3cde0",
showCurve = FALSE,
cex = 1,
main = levels(temp.data$condition)[i]
)
}### High vs. low variability ###
exp1.Qdisgust.posterior.test.var <- filter(
exp1.Qdisgust.posterior.test,
condition == c(
"diff.anger.highVSlow", "diff.neutral.highVSlow", "diff.happiness.highVSlow", "diff.disgust.highVSlow"
)
) %>%
droplevels()
par(mfrow = c(2, 2))
for(i in 1:length(levels(exp1.Qdisgust.posterior.test.var$condition))) {
temp.data <- filter(exp1.Qdisgust.posterior.test.var, condition == levels(condition)[i])
plotPost(
pull(temp.data, ratings),
credMass = 0.95,
compVal = 0,
ROPE = NULL,
xlab = "",
col = "#b3cde0",
showCurve = FALSE,
cex = 1,
main = levels(temp.data$condition)[i]
)
}# subset data
exp2.Qangry.ratings <- ratings %>%
filter(exp == "exp2", question == "angry") %>%
droplevels() %>%
mutate(emotion = factor(emotion,
levels = c("anger", "neutral", "happiness", "disgust") # sort anger first
))
exp2.Qangry.model <- brm(ratings ~ emotion : variability
+ (condition | participant),
data = exp2.Qangry.ratings,
family = cumulative("probit"),
prior = weak.priors,
sample_prior = TRUE,
inits = "0",
control = list(
adapt_delta = .99,
max_treedepth = 15
),
chains = num.chains,
iter = num.iter,
warmup = num.warmup,
thin = num.thin,
algorithm = "sampling",
cores = num.chains,
seed = seed.smorfia
)
saveRDS(exp2.Qangry.model, file = "exp2_Qangry_model.rds", compress = "gzip")# subset data
exp2.Qangry.ratings <- ratings %>%
filter(exp == "exp2", question == "angry") %>%
droplevels() %>%
mutate(emotion = factor(emotion,
levels = c("anger", "neutral", "happiness", "disgust") # sort anger first
))
# load model
exp2.Qangry.model <- readRDS("exp2_Qangry_model.rds")
# summary posterior distributions & diagnostics
summary.exp2.Qangry.model <-
describe_posterior(exp2.Qangry.model,
centrality = "mean",
ci = 0.95,
ci_method = "hdi",
test = NULL,
diagnostic = "all"
) %>%
as_tibble() %>%
separate(Parameter, c("coef", "Parameter"), sep = "_") %>%
dplyr::select(-coef) %>%
mutate(Parameter = recode(factor(Parameter),
"Intercept.1." = "Intercept[1]",
"Intercept.2." = "Intercept[2]",
"Intercept.3." = "Intercept[3]",
"Intercept.4." = "Intercept[4]",
"Intercept.5." = "Intercept[5]",
"Intercept.6." = "Intercept[6]",
"Intercept.7." = "Intercept[7]",
"Intercept.8." = "Intercept[8]",
"emotionanger.variabilityhigh" = "anger, high variability",
"emotionneutral.variabilityhigh" = "neutral, high variability",
"emotionhappiness.variabilityhigh" = "happiness, high variability",
"emotiondisgust.variabilityhigh" = "disgust, high variability",
"emotionanger.variabilitylow" = "anger, low variability",
"emotionneutral.variabilitylow" = "neutral, low variability",
"emotionhappiness.variabilitylow" = "happiness, low variability",
"emotiondisgust.variabilitylow" = "disgust, low variability"
)) %>%
dplyr::select(Parameter, Mean, HDI95_Low = CI_low, HDI95_High = CI_high, Estimated_Sample_Size = ESS, MonteCarlo_StdErr = MCSE)
kable(summary.exp2.Qangry.model, digits = 2)| Parameter | Mean | HDI95_Low | HDI95_High | Estimated_Sample_Size | MonteCarlo_StdErr |
|---|---|---|---|---|---|
| Intercept[1] | -0.85 | -2.94 | 1.34 | 6243 | 0.01 |
| Intercept[2] | -0.12 | -2.26 | 2.02 | 6362 | 0.01 |
| Intercept[3] | 0.44 | -1.66 | 2.62 | 6361 | 0.01 |
| Intercept[4] | 0.90 | -1.20 | 3.08 | 6286 | 0.01 |
| Intercept[5] | 1.77 | -0.33 | 4.01 | 6081 | 0.01 |
| Intercept[6] | 2.45 | 0.25 | 4.62 | 5836 | 0.01 |
| Intercept[7] | 3.32 | 1.07 | 5.53 | 5711 | 0.02 |
| Intercept[8] | 4.28 | 2.09 | 6.70 | 5549 | 0.02 |
| anger, high variability | 4.56 | 2.25 | 6.84 | 5539 | 0.02 |
| neutral, high variability | 0.69 | -1.50 | 2.83 | 6383 | 0.01 |
| happiness, high variability | -3.70 | -6.66 | -0.79 | 9174 | 0.02 |
| disgust, high variability | 0.42 | -1.90 | 2.51 | 6832 | 0.01 |
| anger, low variability | 2.39 | 0.21 | 4.56 | 5951 | 0.01 |
| neutral, low variability | -0.36 | -2.40 | 1.89 | 6419 | 0.01 |
| happiness, low variability | -4.59 | -7.96 | -1.48 | 8720 | 0.02 |
| disgust, low variability | 0.43 | -1.76 | 2.53 | 6519 | 0.01 |
# MCMC chains
exp2.Qangry.model %>%
spread_draws(`b_.*`, regex = TRUE) %>%
dplyr::select("Chain" = ".chain",
"Intercept[1]" = "b_Intercept[1]",
"Intercept[2]" = "b_Intercept[2]",
"Intercept[3]" = "b_Intercept[3]",
"Intercept[4]" = "b_Intercept[4]",
"Intercept[5]" = "b_Intercept[5]",
"Intercept[6]" = "b_Intercept[6]",
"Intercept[7]" = "b_Intercept[7]",
"Intercept[8]" = "b_Intercept[8]",
"anger, high variability" = "b_emotionanger:variabilityhigh",
"neutral, high variability" = "b_emotionneutral:variabilityhigh",
"happiness, high variability" = "b_emotionhappiness:variabilityhigh",
"disgust, high variability" = "b_emotiondisgust:variabilityhigh",
"anger, low variability" = "b_emotionanger:variabilitylow",
"neutral, low variability" = "b_emotionneutral:variabilitylow",
"happiness, low variability" = "b_emotionhappiness:variabilitylow",
"disgust, low variability" = "b_emotiondisgust:variabilitylow"
) %>%
as.data.frame() %>%
mcmc_trace(facet_args = list(ncol = 2)) +
geom_vline(xintercept = seq(1, 9, 1),
linetype = "dotted",
colour = "#999999",
size = .8,
alpha = .5
) +
ggtitle("Chain convergence") +
theme_EmoSSR +
theme(legend.position = "none")### observed vs. predicted data ###
exp2.Qangry.ratings %>%
dplyr::select(exp, participant, question, emotion, variability, condition, ratings) %>% # reorder columns
mutate(predicted = posterior_predict(exp2.Qangry.model)[1, ]) %>% # take only first row of sampled ratings (as example)
dplyr::rename(observed = ratings) %>%
gather(key = data_type, value = ratings, observed:predicted) %>%
ggplot(aes(x = condition,
y = ratings,
color = condition,
fill = condition)) +
geom_pirate(
bars = FALSE,
cis = TRUE,
lines = TRUE, lines_params = list(color = "black"),
points = TRUE, points_params = list(shape = 16, size = 5, alpha = .2),
violins = TRUE, violins_params = list(size = 1),
show.legend = FALSE
) +
scale_fill_viridis_d(option = "cividis") +
scale_color_viridis_d(option = "cividis") +
scale_y_continuous(limits = c(1, 9),
breaks = seq(1, 9, 1)) +
geom_hline(yintercept = seq(1, 9, 1),
linetype = "dotted",
colour = "#999999",
size = .8,
alpha = .5) +
facet_wrap( ~ data_type, scales = "free") +
xlab("") +
ggtitle("pilot 2, \"anger\" question ") +
theme_EmoSSRexp2.Qangry.posterior.test <- exp2.Qangry.model %>%
spread_draws(`b_.*`, regex = TRUE) %>%
dplyr::select(
"anger, high variability" = "b_emotionanger:variabilityhigh",
"neutral, high variability" = "b_emotionneutral:variabilityhigh",
"happiness, high variability" = "b_emotionhappiness:variabilityhigh",
"disgust, high variability" = "b_emotiondisgust:variabilityhigh",
"anger, low variability" = "b_emotionanger:variabilitylow",
"neutral, low variability" = "b_emotionneutral:variabilitylow",
"happiness, low variability" = "b_emotionhappiness:variabilitylow",
"disgust, low variability" = "b_emotiondisgust:variabilitylow"
) %>%
mutate(
### emotion (anger minus others)
## high variability
# anger minus neutral
diff.high.angerVSneutral = `anger, high variability` - `neutral, high variability`,
# anger minus happiness
diff.high.angerVShappiness = `anger, high variability` - `happiness, high variability`,
# anger minus disgust
diff.high.angerVSdisgust = `anger, high variability` - `disgust, high variability`,
## low variability
# anger minus neutral
diff.low.angerVSneutral = `anger, low variability` - `neutral, low variability`,
# anger minus happiness
diff.low.angerVShappiness = `anger, low variability` - `happiness, low variability`,
# anger minus disgust
diff.low.angerVSdisgust = `anger, low variability` - `disgust, low variability`,
### variability (high minus low)
# anger
diff.anger.highVSlow = `anger, high variability` - `anger, low variability`,
# neutral
diff.neutral.highVSlow = `neutral, high variability` - `neutral, low variability`,
# happiness
diff.happiness.highVSlow = `happiness, high variability` - `happiness, low variability`,
# disgust
diff.disgust.highVSlow = `disgust, high variability` - `disgust, low variability`
) %>%
gather(key = "condition", value = "ratings") %>%
filter(condition == c(
"diff.high.angerVSneutral",
"diff.high.angerVShappiness",
"diff.high.angerVSdisgust",
"diff.low.angerVSneutral",
"diff.low.angerVShappiness",
"diff.low.angerVSdisgust",
"diff.anger.highVSlow",
"diff.neutral.highVSlow",
"diff.happiness.highVSlow",
"diff.disgust.highVSlow"
)) %>%
mutate(condition = factor(condition,
levels = c(
"diff.high.angerVSneutral", "diff.high.angerVShappiness", "diff.high.angerVSdisgust",
"diff.low.angerVSneutral", "diff.low.angerVShappiness", "diff.low.angerVSdisgust",
"diff.anger.highVSlow", "diff.neutral.highVSlow", "diff.happiness.highVSlow", "diff.disgust.highVSlow"
)
))
# summary
summary.exp2.Qangry.posterior.test <- exp2.Qangry.posterior.test %>%
group_by(condition) %>%
dplyr::summarize(
mean = mean(ratings),
hdi.95.low = hdi(ratings)[1],
hdi.95.high = hdi(ratings)[2],
evidence.ratio = ifelse(mean < 0,
length(which(ratings < 0)) / length(which(ratings > 0)),
length(which(ratings > 0)) / length(which(ratings < 0))
)
) %>%
ungroup()
kable(summary.exp2.Qangry.posterior.test, digits = 2)| condition | mean | hdi.95.low | hdi.95.high | evidence.ratio |
|---|---|---|---|---|
| diff.high.angerVSneutral | 3.87 | 2.88 | 4.96 | Inf |
| diff.high.angerVShappiness | 8.22 | 6.03 | 11.20 | Inf |
| diff.high.angerVSdisgust | 4.15 | 2.86 | 5.44 | Inf |
| diff.low.angerVSneutral | 2.77 | 1.88 | 3.63 | Inf |
| diff.low.angerVShappiness | 6.94 | 4.19 | 9.75 | Inf |
| diff.low.angerVSdisgust | 1.96 | 1.21 | 2.85 | Inf |
| diff.anger.highVSlow | 2.17 | 1.47 | 2.99 | Inf |
| diff.neutral.highVSlow | 1.04 | 0.25 | 1.76 | 319.00 |
| diff.happiness.highVSlow | 0.87 | -2.62 | 4.57 | 2.08 |
| diff.disgust.highVSlow | -0.01 | -0.96 | 0.95 | 1.04 |
### Anger vs. others ###
exp2.Qangry.posterior.test.emoref <- filter(
exp2.Qangry.posterior.test,
condition == c(
"diff.high.angerVSneutral", "diff.high.angerVShappiness", "diff.high.angerVSdisgust",
"diff.low.angerVSneutral", "diff.low.angerVShappiness", "diff.low.angerVSdisgust"
)
) %>%
droplevels()
par(mfrow = c(3, 2))
for(i in 1:length(levels(exp2.Qangry.posterior.test.emoref$condition))) {
temp.data <- filter(exp2.Qangry.posterior.test.emoref, condition == levels(condition)[i])
plotPost(
pull(temp.data, ratings),
credMass = 0.95,
compVal = 0,
ROPE = NULL,
xlab = "",
col = "#b3cde0",
showCurve = FALSE,
cex = 1,
main = levels(temp.data$condition)[i]
)
}### High vs. low variability ###
exp2.Qangry.posterior.test.var <- filter(
exp2.Qangry.posterior.test,
condition == c(
"diff.anger.highVSlow", "diff.neutral.highVSlow", "diff.happiness.highVSlow", "diff.disgust.highVSlow"
)
) %>%
droplevels()
par(mfrow = c(2, 2))
for(i in 1:length(levels(exp2.Qangry.posterior.test.var$condition))) {
temp.data <- filter(exp2.Qangry.posterior.test.var, condition == levels(condition)[i])
plotPost(
pull(temp.data, ratings),
credMass = 0.95,
compVal = 0,
ROPE = NULL,
xlab = "",
col = "#b3cde0",
showCurve = FALSE,
cex = 1,
main = levels(temp.data$condition)[i]
)
}# subset data
exp2.Qhappy.ratings <- ratings %>%
filter(exp == "exp2", question == "happy") %>%
droplevels() %>%
mutate(emotion = factor(emotion,
levels = c("happiness", "neutral", "anger", "disgust") # sort happiness first
))
exp2.Qhappy.model <- brm(ratings ~ emotion : variability
+ (condition | participant),
data = exp2.Qhappy.ratings,
family = cumulative("probit"),
prior = weak.priors,
sample_prior = TRUE,
inits = "0",
control = list(
adapt_delta = .99,
max_treedepth = 15
),
chains = num.chains,
iter = num.iter,
warmup = num.warmup,
thin = num.thin,
algorithm = "sampling",
cores = num.chains,
seed = seed.smorfia
)
saveRDS(exp2.Qhappy.model, file = "exp2_Qhappy_model.rds", compress = "gzip")# subset data
exp2.Qhappy.ratings <- ratings %>%
filter(exp == "exp2", question == "happy") %>%
droplevels() %>%
mutate(emotion = factor(emotion,
levels = c("happiness", "neutral", "anger", "disgust") # sort happiness first
))
# load model
exp2.Qhappy.model <- readRDS("exp2_Qhappy_model.rds")
# summary posterior distributions & diagnostics
summary.exp2.Qhappy.model <-
describe_posterior(exp2.Qhappy.model,
centrality = "mean",
ci = 0.95,
ci_method = "hdi",
test = NULL,
diagnostic = "all"
) %>%
as_tibble() %>%
separate(Parameter, c("coef", "Parameter"), sep = "_") %>%
dplyr::select(-coef) %>%
mutate(Parameter = recode(factor(Parameter),
"Intercept.1." = "Intercept[1]",
"Intercept.2." = "Intercept[2]",
"Intercept.3." = "Intercept[3]",
"Intercept.4." = "Intercept[4]",
"Intercept.5." = "Intercept[5]",
"Intercept.6." = "Intercept[6]",
"Intercept.7." = "Intercept[7]",
"Intercept.8." = "Intercept[8]",
"emotionanger.variabilityhigh" = "anger, high variability",
"emotionneutral.variabilityhigh" = "neutral, high variability",
"emotionhappiness.variabilityhigh" = "happiness, high variability",
"emotiondisgust.variabilityhigh" = "disgust, high variability",
"emotionanger.variabilitylow" = "anger, low variability",
"emotionneutral.variabilitylow" = "neutral, low variability",
"emotionhappiness.variabilitylow" = "happiness, low variability",
"emotiondisgust.variabilitylow" = "disgust, low variability"
)) %>%
dplyr::select(Parameter, Mean, HDI95_Low = CI_low, HDI95_High = CI_high, Estimated_Sample_Size = ESS, MonteCarlo_StdErr = MCSE)
kable(summary.exp2.Qhappy.model, digits = 2)| Parameter | Mean | HDI95_Low | HDI95_High | Estimated_Sample_Size | MonteCarlo_StdErr |
|---|---|---|---|---|---|
| Intercept[1] | -1.51 | -3.73 | 0.65 | 5492 | 0.02 |
| Intercept[2] | -0.57 | -2.74 | 1.62 | 5637 | 0.01 |
| Intercept[3] | 0.26 | -1.91 | 2.43 | 5720 | 0.01 |
| Intercept[4] | 0.65 | -1.51 | 2.86 | 5879 | 0.01 |
| Intercept[5] | 2.12 | -0.14 | 4.30 | 5996 | 0.01 |
| Intercept[6] | 2.74 | 0.52 | 5.04 | 6054 | 0.01 |
| Intercept[7] | 3.95 | 1.58 | 6.32 | 5886 | 0.02 |
| Intercept[8] | 5.85 | 3.26 | 8.57 | 5132 | 0.02 |
| happiness, high variability | 5.55 | 2.98 | 8.19 | 5577 | 0.02 |
| neutral, high variability | -0.67 | -2.87 | 1.65 | 5954 | 0.02 |
| anger, high variability | -2.71 | -5.01 | -0.57 | 5506 | 0.02 |
| disgust, high variability | -3.02 | -5.26 | -0.66 | 5617 | 0.02 |
| happiness, low variability | 5.08 | 2.46 | 7.62 | 5763 | 0.02 |
| neutral, low variability | 0.19 | -1.96 | 2.36 | 5764 | 0.01 |
| anger, low variability | -1.78 | -4.00 | 0.41 | 5573 | 0.02 |
| disgust, low variability | -2.75 | -5.01 | -0.47 | 5530 | 0.02 |
# MCMC chains
exp2.Qhappy.model %>%
spread_draws(`b_.*`, regex = TRUE) %>%
dplyr::select("Chain" = ".chain",
"Intercept[1]" = "b_Intercept[1]",
"Intercept[2]" = "b_Intercept[2]",
"Intercept[3]" = "b_Intercept[3]",
"Intercept[4]" = "b_Intercept[4]",
"Intercept[5]" = "b_Intercept[5]",
"Intercept[6]" = "b_Intercept[6]",
"Intercept[7]" = "b_Intercept[7]",
"Intercept[8]" = "b_Intercept[8]",
"anger, high variability" = "b_emotionanger:variabilityhigh",
"neutral, high variability" = "b_emotionneutral:variabilityhigh",
"happiness, high variability" = "b_emotionhappiness:variabilityhigh",
"disgust, high variability" = "b_emotiondisgust:variabilityhigh",
"anger, low variability" = "b_emotionanger:variabilitylow",
"neutral, low variability" = "b_emotionneutral:variabilitylow",
"happiness, low variability" = "b_emotionhappiness:variabilitylow",
"disgust, low variability" = "b_emotiondisgust:variabilitylow"
) %>%
as.data.frame() %>%
mcmc_trace(facet_args = list(ncol = 2)) +
geom_vline(xintercept = seq(1, 9, 1),
linetype = "dotted",
colour = "#999999",
size = .8,
alpha = .5
) +
ggtitle("Chain convergence") +
theme_EmoSSR +
theme(legend.position = "none")### observed vs. predicted data ###
exp2.Qhappy.ratings %>%
dplyr::select(exp, participant, question, emotion, variability, condition, ratings) %>% # reorder columns
mutate(predicted = posterior_predict(exp2.Qhappy.model)[1, ]) %>% # take only first row of sampled ratings (as example)
dplyr::rename(observed = ratings) %>%
gather(key = data_type, value = ratings, observed:predicted) %>%
ggplot(aes(x = condition,
y = ratings,
color = condition,
fill = condition)) +
geom_pirate(
bars = FALSE,
cis = TRUE,
lines = TRUE, lines_params = list(color = "black"),
points = TRUE, points_params = list(shape = 16, size = 5, alpha = .2),
violins = TRUE, violins_params = list(size = 1),
show.legend = FALSE
) +
scale_fill_viridis_d(option = "cividis") +
scale_color_viridis_d(option = "cividis") +
scale_y_continuous(limits = c(1, 9),
breaks = seq(1, 9, 1)) +
geom_hline(yintercept = seq(1, 9, 1),
linetype = "dotted",
colour = "#999999",
size = .8,
alpha = .5) +
facet_wrap( ~ data_type, scales = "free") +
xlab("") +
ggtitle("experiment 2, \"happiness\" question ") +
theme_EmoSSRexp2.Qhappy.posterior.test <- exp2.Qhappy.model %>%
spread_draws(`b_.*`, regex = TRUE) %>%
dplyr::select(
"anger, high variability" = "b_emotionanger:variabilityhigh",
"neutral, high variability" = "b_emotionneutral:variabilityhigh",
"happiness, high variability" = "b_emotionhappiness:variabilityhigh",
"disgust, high variability" = "b_emotiondisgust:variabilityhigh",
"anger, low variability" = "b_emotionanger:variabilitylow",
"neutral, low variability" = "b_emotionneutral:variabilitylow",
"happiness, low variability" = "b_emotionhappiness:variabilitylow",
"disgust, low variability" = "b_emotiondisgust:variabilitylow"
) %>%
mutate(
### emotion (happiness minus others)
## high variability
# happiness minus neutral
diff.high.happinessVSneutral = `happiness, high variability` - `neutral, high variability`,
# happiness minus anger
diff.high.happinessVSanger = `happiness, high variability` - `anger, high variability`,
# happiness minus disgust
diff.high.happinessVSdisgust = `happiness, high variability` - `disgust, high variability`,
## low variability
# happiness minus neutral
diff.low.happinessVSneutral = `happiness, low variability` - `neutral, low variability`,
# happiness minus anger
diff.low.happinessVSanger = `happiness, low variability` - `anger, low variability`,
# happiness minus disgust
diff.low.happinessVSdisgust = `happiness, low variability` - `disgust, low variability`,
### variability (high minus low)
# anger
diff.anger.highVSlow = `anger, high variability` - `anger, low variability`,
# neutral
diff.neutral.highVSlow = `neutral, high variability` - `neutral, low variability`,
# happiness
diff.happiness.highVSlow = `happiness, high variability` - `happiness, low variability`,
# disgust
diff.disgust.highVSlow = `disgust, high variability` - `disgust, low variability`
) %>%
gather(key = "condition", value = "ratings") %>%
filter(condition == c(
"diff.high.happinessVSneutral",
"diff.high.happinessVSanger",
"diff.high.happinessVSdisgust",
"diff.low.happinessVSneutral",
"diff.low.happinessVSanger",
"diff.low.happinessVSdisgust",
"diff.anger.highVSlow",
"diff.neutral.highVSlow",
"diff.happiness.highVSlow",
"diff.disgust.highVSlow"
)) %>%
mutate(condition = factor(condition,
levels = c(
"diff.high.happinessVSneutral", "diff.high.happinessVSanger", "diff.high.happinessVSdisgust",
"diff.low.happinessVSneutral", "diff.low.happinessVSanger", "diff.low.happinessVSdisgust",
"diff.anger.highVSlow", "diff.neutral.highVSlow", "diff.happiness.highVSlow", "diff.disgust.highVSlow"
)
))
# summary
summary.exp2.Qhappy.posterior.test <- exp2.Qhappy.posterior.test %>%
group_by(condition) %>%
dplyr::summarize(
mean = mean(ratings),
hdi.95.low = hdi(ratings)[1],
hdi.95.high = hdi(ratings)[2],
evidence.ratio = ifelse(mean < 0,
length(which(ratings < 0)) / length(which(ratings > 0)),
length(which(ratings > 0)) / length(which(ratings < 0))
)
) %>%
ungroup()
kable(summary.exp2.Qhappy.posterior.test, digits = 2)| condition | mean | hdi.95.low | hdi.95.high | evidence.ratio |
|---|---|---|---|---|
| diff.high.happinessVSneutral | 6.23 | 4.18 | 8.51 | Inf |
| diff.high.happinessVSanger | 8.25 | 6.11 | 10.62 | Inf |
| diff.high.happinessVSdisgust | 8.55 | 6.35 | 11.02 | Inf |
| diff.low.happinessVSneutral | 4.90 | 3.09 | 6.75 | Inf |
| diff.low.happinessVSanger | 6.88 | 4.81 | 8.79 | Inf |
| diff.low.happinessVSdisgust | 7.81 | 5.59 | 9.87 | Inf |
| diff.anger.highVSlow | -0.93 | -1.70 | -0.07 | 71.73 |
| diff.neutral.highVSlow | -0.85 | -1.74 | 0.29 | 24.00 |
| diff.happiness.highVSlow | 0.46 | -0.99 | 1.86 | 3.07 |
| diff.disgust.highVSlow | -0.27 | -1.55 | 0.77 | 2.26 |
### Happiness vs. others ###
exp2.Qhappy.posterior.test.emoref <- filter(
exp2.Qhappy.posterior.test,
condition == c(
"diff.high.happinessVSneutral", "diff.high.happinessVSanger", "diff.high.happinessVSdisgust",
"diff.low.happinessVSneutral", "diff.low.happinessVSanger", "diff.low.happinessVSdisgust"
)
) %>%
droplevels()
par(mfrow = c(3, 2))
for(i in 1:length(levels(exp2.Qhappy.posterior.test.emoref$condition))) {
temp.data <- filter(exp2.Qhappy.posterior.test.emoref, condition == levels(condition)[i])
plotPost(
pull(temp.data, ratings),
credMass = 0.95,
compVal = 0,
ROPE = NULL,
xlab = "",
col = "#b3cde0",
showCurve = FALSE,
cex = 1,
main = levels(temp.data$condition)[i]
)
}### High vs. low variability ###
exp2.Qhappy.posterior.test.var <- filter(
exp2.Qhappy.posterior.test,
condition == c(
"diff.anger.highVSlow", "diff.neutral.highVSlow", "diff.happiness.highVSlow", "diff.disgust.highVSlow"
)
) %>%
droplevels()
par(mfrow = c(2, 2))
for(i in 1:length(levels(exp2.Qhappy.posterior.test.var$condition))) {
temp.data <- filter(exp2.Qhappy.posterior.test.var, condition == levels(condition)[i])
plotPost(
pull(temp.data, ratings),
credMass = 0.95,
compVal = 0,
ROPE = NULL,
xlab = "",
col = "#b3cde0",
showCurve = FALSE,
cex = 1,
main = levels(temp.data$condition)[i]
)
}# subset data
exp2.Qdisgust.ratings <- ratings %>%
filter(exp == "exp2", question == "disgusted") %>%
droplevels() %>%
mutate(emotion = factor(emotion,
levels = c("disgust", "neutral", "anger", "happiness") # sort disgust first
))
exp2.Qdisgust.model <- brm(ratings ~ emotion : variability
+ (condition | participant),
data = exp2.Qdisgust.ratings,
family = cumulative("probit"),
prior = weak.priors,
sample_prior = TRUE,
inits = "0",
control = list(
adapt_delta = .99,
max_treedepth = 15
),
chains = num.chains,
iter = num.iter,
warmup = num.warmup,
thin = num.thin,
algorithm = "sampling",
cores = num.chains,
seed = seed.smorfia
)
saveRDS(exp2.Qdisgust.model, file = "exp2_Qdisgust_model.rds", compress = "gzip")# subset data
exp2.Qdisgust.ratings <- ratings %>%
filter(exp == "exp2", question == "disgusted") %>%
droplevels() %>%
mutate(emotion = factor(emotion,
levels = c("disgust", "neutral", "anger", "happiness") # sort disgust first
))
# load model
exp2.Qdisgust.model <- readRDS("exp2_Qdisgust_model.rds")
# summary posterior distributions & diagnostics
summary.exp2.Qdisgust.model <-
describe_posterior(exp2.Qdisgust.model,
centrality = "mean",
ci = 0.95,
ci_method = "hdi",
test = NULL,
diagnostic = "all"
) %>%
as_tibble() %>%
separate(Parameter, c("coef", "Parameter"), sep = "_") %>%
dplyr::select(-coef) %>%
mutate(Parameter = recode(factor(Parameter),
"Intercept.1." = "Intercept[1]",
"Intercept.2." = "Intercept[2]",
"Intercept.3." = "Intercept[3]",
"Intercept.4." = "Intercept[4]",
"Intercept.5." = "Intercept[5]",
"Intercept.6." = "Intercept[6]",
"Intercept.7." = "Intercept[7]",
"Intercept.8." = "Intercept[8]",
"emotionanger.variabilityhigh" = "anger, high variability",
"emotionneutral.variabilityhigh" = "neutral, high variability",
"emotionhappiness.variabilityhigh" = "happiness, high variability",
"emotiondisgust.variabilityhigh" = "disgust, high variability",
"emotionanger.variabilitylow" = "anger, low variability",
"emotionneutral.variabilitylow" = "neutral, low variability",
"emotionhappiness.variabilitylow" = "happiness, low variability",
"emotiondisgust.variabilitylow" = "disgust, low variability"
)) %>%
dplyr::select(Parameter, Mean, HDI95_Low = CI_low, HDI95_High = CI_high, Estimated_Sample_Size = ESS, MonteCarlo_StdErr = MCSE)
kable(summary.exp2.Qdisgust.model, digits = 2)| Parameter | Mean | HDI95_Low | HDI95_High | Estimated_Sample_Size | MonteCarlo_StdErr |
|---|---|---|---|---|---|
| Intercept[1] | -1.44 | -3.55 | 0.72 | 7608 | 0.01 |
| Intercept[2] | -0.39 | -2.48 | 1.78 | 7759 | 0.01 |
| Intercept[3] | 0.15 | -2.04 | 2.24 | 7783 | 0.01 |
| Intercept[4] | 0.46 | -1.59 | 2.68 | 7819 | 0.01 |
| Intercept[5] | 1.18 | -0.95 | 3.34 | 7805 | 0.01 |
| Intercept[6] | 1.71 | -0.35 | 3.98 | 7765 | 0.01 |
| Intercept[7] | 2.23 | 0.00 | 4.36 | 7664 | 0.01 |
| Intercept[8] | 3.51 | 1.15 | 5.71 | 6736 | 0.01 |
| disgust, high variability | 3.99 | 1.54 | 6.51 | 7636 | 0.01 |
| neutral, high variability | -0.05 | -2.18 | 2.16 | 8025 | 0.01 |
| anger, high variability | -0.48 | -2.67 | 1.60 | 7827 | 0.01 |
| happiness, high variability | -3.61 | -6.31 | -0.79 | 9212 | 0.01 |
| disgust, low variability | 3.07 | 0.73 | 5.31 | 7593 | 0.01 |
| neutral, low variability | -1.11 | -3.27 | 1.06 | 7765 | 0.01 |
| anger, low variability | 0.49 | -1.66 | 2.60 | 7876 | 0.01 |
| happiness, low variability | -2.29 | -4.48 | -0.01 | 7681 | 0.01 |
# MCMC chains
exp2.Qdisgust.model %>%
spread_draws(`b_.*`, regex = TRUE) %>%
dplyr::select("Chain" = ".chain",
"Intercept[1]" = "b_Intercept[1]",
"Intercept[2]" = "b_Intercept[2]",
"Intercept[3]" = "b_Intercept[3]",
"Intercept[4]" = "b_Intercept[4]",
"Intercept[5]" = "b_Intercept[5]",
"Intercept[6]" = "b_Intercept[6]",
"Intercept[7]" = "b_Intercept[7]",
"Intercept[8]" = "b_Intercept[8]",
"anger, high variability" = "b_emotionanger:variabilityhigh",
"neutral, high variability" = "b_emotionneutral:variabilityhigh",
"happiness, high variability" = "b_emotionhappiness:variabilityhigh",
"disgust, high variability" = "b_emotiondisgust:variabilityhigh",
"anger, low variability" = "b_emotionanger:variabilitylow",
"neutral, low variability" = "b_emotionneutral:variabilitylow",
"happiness, low variability" = "b_emotionhappiness:variabilitylow",
"disgust, low variability" = "b_emotiondisgust:variabilitylow"
) %>%
as.data.frame() %>%
mcmc_trace(facet_args = list(ncol = 2)) +
geom_vline(xintercept = seq(1, 9, 1),
linetype = "dotted",
colour = "#999999",
size = .8,
alpha = .5
) +
ggtitle("Chain convergence") +
theme_EmoSSR +
theme(legend.position = "none")### observed vs. predicted data ###
exp2.Qdisgust.ratings %>%
dplyr::select(exp, participant, question, emotion, variability, condition, ratings) %>% # reorder columns
mutate(predicted = posterior_predict(exp2.Qdisgust.model)[1, ]) %>% # take only first row of sampled ratings (as example)
dplyr::rename(observed = ratings) %>%
gather(key = data_type, value = ratings, observed:predicted) %>%
ggplot(aes(x = condition,
y = ratings,
color = condition,
fill = condition)) +
geom_pirate(
bars = FALSE,
cis = TRUE,
lines = TRUE, lines_params = list(color = "black"),
points = TRUE, points_params = list(shape = 16, size = 5, alpha = .2),
violins = TRUE, violins_params = list(size = 1),
show.legend = FALSE
) +
scale_fill_viridis_d(option = "cividis") +
scale_color_viridis_d(option = "cividis") +
scale_y_continuous(limits = c(1, 9),
breaks = seq(1, 9, 1)) +
geom_hline(yintercept = seq(1, 9, 1),
linetype = "dotted",
colour = "#999999",
size = .8,
alpha = .5) +
facet_wrap( ~ data_type, scales = "free") +
xlab("") +
ggtitle("experiment 2, \"disgusted\" question ") +
theme_EmoSSRexp2.Qdisgust.posterior.test <- exp2.Qdisgust.model %>%
spread_draws(`b_.*`, regex = TRUE) %>%
dplyr::select(
"anger, high variability" = "b_emotionanger:variabilityhigh",
"neutral, high variability" = "b_emotionneutral:variabilityhigh",
"happiness, high variability" = "b_emotionhappiness:variabilityhigh",
"disgust, high variability" = "b_emotiondisgust:variabilityhigh",
"anger, low variability" = "b_emotionanger:variabilitylow",
"neutral, low variability" = "b_emotionneutral:variabilitylow",
"happiness, low variability" = "b_emotionhappiness:variabilitylow",
"disgust, low variability" = "b_emotiondisgust:variabilitylow"
) %>%
mutate(
### emotion (disgust minus others)
## high variability
# disgust minus neutral
diff.high.disgustVSneutral = `disgust, high variability` - `neutral, high variability`,
# disgust minus anger
diff.high.disgustVSanger = `disgust, high variability` - `anger, high variability`,
# disgust minus happiness
diff.high.disgustVShappiness = `disgust, high variability` - `happiness, high variability`,
## low variability
# disgust minus neutral
diff.low.disgustVSneutral = `disgust, low variability` - `neutral, low variability`,
# disgust minus anger
diff.low.disgustVSanger = `disgust, low variability` - `anger, low variability`,
# disgust minus happiness
diff.low.disgustVShappiness = `disgust, low variability` - `happiness, low variability`,
### variability (high minus low)
# anger
diff.anger.highVSlow = `anger, high variability` - `anger, low variability`,
# neutral
diff.neutral.highVSlow = `neutral, high variability` - `neutral, low variability`,
# happiness
diff.happiness.highVSlow = `happiness, high variability` - `happiness, low variability`,
# disgust
diff.disgust.highVSlow = `disgust, high variability` - `disgust, low variability`
) %>%
gather(key = "condition", value = "ratings") %>%
filter(condition == c(
"diff.high.disgustVSneutral",
"diff.high.disgustVSanger",
"diff.high.disgustVShappiness",
"diff.low.disgustVSneutral",
"diff.low.disgustVSanger",
"diff.low.disgustVShappiness",
"diff.anger.highVSlow",
"diff.neutral.highVSlow",
"diff.happiness.highVSlow",
"diff.disgust.highVSlow"
)) %>%
mutate(condition = factor(condition,
levels = c(
"diff.high.disgustVSneutral", "diff.high.disgustVSanger", "diff.high.disgustVShappiness",
"diff.low.disgustVSneutral", "diff.low.disgustVSanger", "diff.low.disgustVShappiness",
"diff.anger.highVSlow", "diff.neutral.highVSlow", "diff.happiness.highVSlow", "diff.disgust.highVSlow"
)
))
# summary
summary.exp2.Qdisgust.posterior.test <- exp2.Qdisgust.posterior.test %>%
group_by(condition) %>%
dplyr::summarize(
mean = mean(ratings),
hdi.95.low = hdi(ratings)[1],
hdi.95.high = hdi(ratings)[2],
evidence.ratio = ifelse(mean < 0,
length(which(ratings < 0)) / length(which(ratings > 0)),
length(which(ratings > 0)) / length(which(ratings < 0))
)
) %>%
ungroup()
kable(summary.exp2.Qdisgust.posterior.test, digits = 2)| condition | mean | hdi.95.low | hdi.95.high | evidence.ratio |
|---|---|---|---|---|
| diff.high.disgustVSneutral | 4.02 | 2.32 | 5.57 | Inf |
| diff.high.disgustVSanger | 4.49 | 3.06 | 6.19 | Inf |
| diff.high.disgustVShappiness | 7.60 | 5.05 | 10.17 | Inf |
| diff.low.disgustVSneutral | 4.20 | 2.93 | 5.57 | Inf |
| diff.low.disgustVSanger | 2.60 | 1.52 | 3.60 | Inf |
| diff.low.disgustVShappiness | 5.38 | 3.92 | 6.94 | Inf |
| diff.anger.highVSlow | -0.97 | -1.56 | -0.41 | 799.00 |
| diff.neutral.highVSlow | 1.07 | 0.19 | 1.86 | 159.00 |
| diff.happiness.highVSlow | -1.30 | -3.25 | 0.66 | 10.85 |
| diff.disgust.highVSlow | 0.93 | -0.30 | 2.33 | 16.20 |
### Disgust vs. others ###
exp2.Qdisgust.posterior.test.emoref <- filter(
exp2.Qdisgust.posterior.test,
condition == c(
"diff.high.disgustVSneutral", "diff.high.disgustVSanger", "diff.high.disgustVShappiness",
"diff.low.disgustVSneutral", "diff.low.disgustVSanger", "diff.low.disgustVShappiness"
)
) %>%
droplevels()
par(mfrow = c(3, 2))
for(i in 1:length(levels(exp2.Qdisgust.posterior.test.emoref$condition))) {
temp.data <- filter(exp2.Qdisgust.posterior.test.emoref, condition == levels(condition)[i])
plotPost(
pull(temp.data, ratings),
credMass = 0.95,
compVal = 0,
ROPE = NULL,
xlab = "",
col = "#b3cde0",
showCurve = FALSE,
cex = 1,
main = levels(temp.data$condition)[i]
)
}### High vs. low variability ###
exp2.Qdisgust.posterior.test.var <- filter(
exp2.Qdisgust.posterior.test,
condition == c(
"diff.anger.highVSlow", "diff.neutral.highVSlow", "diff.happiness.highVSlow", "diff.disgust.highVSlow"
)
) %>%
droplevels()
par(mfrow = c(2, 2))
for(i in 1:length(levels(exp2.Qdisgust.posterior.test.var$condition))) {
temp.data <- filter(exp2.Qdisgust.posterior.test.var, condition == levels(condition)[i])
plotPost(
pull(temp.data, ratings),
credMass = 0.95,
compVal = 0,
ROPE = NULL,
xlab = "",
col = "#b3cde0",
showCurve = FALSE,
cex = 1,
main = levels(temp.data$condition)[i]
)
}